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SUPERSONIC  FLOW  OVER  BODIES  OF  REVOLUTION 
(WITH  SPECIAL  REFERENCE  TO  HIGH  SPEED  COMPUTING) 

ABSTRACT 

With  the  advent  of  large-scale  high-speed  computing  machines,  it  has 
become  feasible  to  solve  certain  supersonic  flow  problems  by  numerical 
methods  using  the  exact  hydrodynamical  equations  instead  of  resorting  to 
linearization  or  graphical  methods.  This  report  describes  in  detail  one 
such  numerical  method;  namely,  an  efficient  form  of  the  method  of  charac¬ 
teristics.--. 

■ — Characteristic  equations  are  derived  for  supersonic,  steady,  invis- 
cid,  isoenergetic  flows  in  terms  of  a  variety  of  dependent  variables .y  The 
computation  described  is  applicable  to  non-yawed  bodies  of  revolutiop/  hav¬ 
ing  pointed  noses  and  fairly  arbitrary  contours,  which  lie  in^^wtfform 
stream  moving  fast  enough  to  produce  a  shock-wave  at  tiie-'-ndSeand  maintain 
supersonic  flow  everywhere.  The  c omputational^.pirrcS^ure  is  divided  into 
several  parts:  Taylor-Maccoll,  cojgjea^ncfitour,  general,  and  shock  proc¬ 
esses.  Equations  and  bo^ndary^Conditions  are  given  for  each  of  these 
procedures 

discussion  is  given  of  several  methods  of  numerically  solving  sys¬ 
tems  of  1st  order  ordinary  differential  equations,  such  as  are  encountered 
in  the  Taylor-Maccoll  and  corner  processes.  The  other  computations  involve 
approximating  partial  derivatives  by  difference  quotients  and  solving  on 
a  finite  grid  of  points.  Solutions  are  derived  for  the  cases  in  which  the 
derivatives  are  approximated  to  1st,  2nd,  and  3rd  orders. 

An  empirical  study  is  made  of  the  error  due  to  the  introduction  of 
finite  differences.  This  is  based  on  the  results  of  a  particular  calcu¬ 
lation  performed  on  the  ENIAC.  It  is  shown  that  a  knowledge  of  the  cap¬ 
ture  of  the  error?  leads  to  a  procedure  for  extrapolation  to  zero  grid 
size,  which  reduces  by  a  factor  of  ten  the  total  labor  required  to  ob¬ 
tain  a  solution  correct  to  about  four  significant  figures. 
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SUPERSONIC  FLOW  OVER  BODIES  OF  REVOLUTION 
(WITH  SPECIAL  REFERENCE  TO  HIGH  SPEED  COMPUTING) 


Section  1.  Introduction 

The  mathematical  basis  for  computing  the  velocity,  density,  and 
pressure  distribution  of  air  flowing  faster  than  sound  over  plane  bodies 

11  12 

and  bodies  of  revolution  has  been  laid  by  Riemann  *  ,  Picard  *  ,  Hada- 

mard^*-^,  Goursat'*'  *^,  Lewy^*^,  Friedrichs  and  Lewy^-*^,  Frankl  and 
17  18 

Aleksieva  *  ,  Courant  and  lax  *  ,  and  others  .  The  problem  is  that  of 
solving  a  non-linear  system  of  hyperbolic  partial  differential  equations 
with  boundary  conditions  given  on  a  known  curve  (the  bod^r  contour)  and 
on  a  curve  not  known  in  advance  (the  shock  wave). 

Methods  have  been  known  for  fifteen  years  for  solving  the  exact 
equations  (without  friction,  with  rotation)  for  supersonic  flow  about 
plane  and  axial  bodies.  Heretofore  only  slight  use  has  been  made  of 
them,  however,  because  of  the  extreme  tediousness  of  the  numerical  com¬ 
putations.  Instead,  the  solution  of  supersonic  flow  problems  has  pro¬ 
ceeded  along  two  main  lines:  (1)  graphical  and  3emtgraphical  procedures, 
developed  especially  by  Prandtl,  Busemann,  Sauer,  Tollmien,  Guderley, 
and  others  of  the  German  school;  and  (2)  linearization  of  the  hydrody- 
namical  equations.  linear  problems  are  easier  to  solve;  whole  families 
of  solutions  may  often  be  obtained  exhibiting  the  variation  of  the  solu¬ 
tions  with  important  parameters.  Indeed,  even  if  the  linear  problem  has 
been  obtained  by  neglecting  some  moderately  large  terms,  the  solution  is 
often  very  valuable  qualitatively  in  guiding  the  intuition. 

1.1  R.  Courant  and  D.  Hilbert,  Methoden  der  Mathematischen  Fhysik,II; 
p.  311*  Julius  Springer,  Berlin,  1937. 

1.2  Es  Picard,  Trait e  d1 analyse,  II.  Paris  (3rd  ed.  1926) 

1.3  J .  Hadamard,  Lecons  sur  le  Problem  de  Cauchy;  p.  487.  Paris,  1932 

1.4  E.  Goursat,  Cours  d*Analyse,  11;  p.  360.  Paris  (4th  ed.  1924) 

1.5  K.  Lewy,  ’’Ueber  das  Anfangswerfc problem  bei  einer  hyperbolischen 

nichtlinearen  partiellen  Differentialgleichung  zweiter  Ordnung  mit 
zwei  unabhangigen  VerA'nderlichen,  "Mathematische  Annalen,  vol.  98 
(1927),  pp.  179-191.  " 

1.6  K.  Friedrichs  and  H.  Lewy,  ’'Das  Anfangswertproblem  einer  beliebigen 
nichtlinearen  hyperbolischen  Differentialgleichung  beliebiger  Ord¬ 
nung  in  zwei  Variabeln.  Existenz,  Eindeutigkeit  und  Abhangigkeits- 
bereich  der  LSsung,’*  Mat hematis che .  Annalen,  vol.  99  (1928),' pp. 
220-221. 

1.7  F.  Frankl  and  P.  Aleksieva,  "Two  Boundary  Value  Problems  from  the 
Theory  of  Hyperbolic  Partial  Differential.  Equations  with  Applica¬ 
tions  to  Supersonic  Gas  Flow",  Rec.  Math.  Mosc.,  T.  41:3  (1934)* 
(Also  BRL  Report  X-123;  Aberdeen  Proving  Ground,  Maryland.) 

1.8  R.  Courant  and  P.  Lax,  rt0n  Nonlinear  Partial  Differential  Equations 
with  Two  Independent  Variables’*,  Communications  on  Pure  and  Applied 
Mathematics,  Vol.  II,  nos.  2-3  (1949);  p. 


5 


PRECEDING  PAGE  BLANK 


With  the  advent  of  high  speed  computing  devices  snch  as  the  ENIAC 
now  operating  at  the  BRL  or  the  EDVAC  being  installed  at  the  BRL,  a  shift 
of  emphasis  will  take  place.  A  greater  effort  will  be  devoted  to  the  solu¬ 
tion  of  the  exact  equations.  It  will  be  possible  to  solve  these  equations 
so  rapidly  that  parameters  may  again  be  introduced.  Since  the  machines 
are  even  able  to  think  in  an  elementary  way,  they  can  be  made  to  solve  in 
a  numerical  manner  such  problems  in  the  Calculus  of  Variations  as  deter¬ 
mining  the  head  shape  of  given  diameter  and  head  length  which  will  lead 
to  minimal  head  drag. 

This  report  has  been  written  in  an  effort  to  accelerate  the  change 
in  emphasis.  It  includes  some  results  obtained  using  the  ENIAC.  It  is 
expected  that  more  ENIAC  and  EDVAC  results  will  appear  in  later  reports. 

Section  2.  Fundamental  Equations 

Introduction 


The  problems  considered  in  this  report  are  all  of  the  following  type: 
air  flows  steadily  and  supersonically,  from  a  region  of  uniformity,  past 
a  body  which  may  be  plane  or  have  symmetry  of  revolution.  If  there  is  a 
3hock  wave,  the  Ifech  number  i3  assumed  large  enough  and  the  initial  flow 
deviation  small  enough  so  that  the  shock  front  is  attached  to  the  body 
at  a  known  point,  and  the  velocity  is  everywhere  supersonic.  Air  is  con¬ 
sidered  a  perfect  gas,  body  forces  and  friction  (therefore  heat  conduction) 
are  neglected,  but  rotation  of  the  flow  caused  by  a  curved  shock  wave  is 
allowed. 


With  these  restrictions  the  continuity,  energy,  and  Euler  equations 

!«1 

(2.1)  V  .  (^  U)  =  0, 

(2.2)  Q2  + - 1 -  A2  =  C2, 

Cr~  i)  - 

(2.3)  (3  .  V )  5  =  -~Vp, 


where  $ » /° ,  p ,Y  »  and  A  are  the  velocity  vector,  density  of  the  air, 
pressure  of  the  air,  ratio  of  specific  heats,  and  velocity  of  sound, 
respectively. 

Equation  (2.2)  shows  that  as  the  velocity  of  sound  approaches  zero 
the  velocity  approaches  a  limit  C.  Equation  (2.2)  holds  across  a  shock 
2  #2 

wave  *  also,  and  therefore  C  is  the  same  for  all  parts  of  the  fluid. 

2.1  R.  Courant  and  R.  0.  Jtiedrichs,  Supersonic  Flow  and  Shock  Ihves; 
pp.  14,  22,  15,  Interscience  Publishers  Inc.,  New  York,  1§48 

2.2  W.  F.  Durand,  Aerodynamic  Theory;  Vol.  m,  p.  217. 
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Accordingly,  we  shall  take  it  as  the  unit  velocity,  setting 


(2.4)  <1  =  Q/C ,  a  “  A/C. 

With  this  notation  equations  (2.1),  (2.2)  and  (2.3)  become 

(2.5)  V  .  (/O  q)  =  0, 

(2.6)  q2  +  a2  =1, 

r- 1 

and  2 

(2.7)  (q  .  7  )  q  =  -  V  F  • 

rp 

For  a  perfect  steady  gas  without  friction  and  body  forces,  it  can 
be  ::hown  that  in  those  regions  where  there  is  no  shock-wave  the  constancy 
of  the  entropy  on  streamlines  follows  from  equations  (2.1),  (2.2)  and 
(2.3).  Hence 

(2.8)  q  .  V  5  =0  where  p/^o  Y  -  e^cv 
or 

(2.9)  p  (q  .  V  p)  =  Y  P  (9  •  V/Q  ) 

Special  Coordinate  Systems 


In  the  case  of  flow  past  or  through  a  body  of  revolution  we  shall 
introduce  the  axis  of  symmetry  a3  the  x-axis  with  the  orientation  of  the 
free  stream  velocity  vector  q^.  We  shall  let  the  y-axis  be  a  line 

through  the  leading  point  A,  perpendicular  to  the  axis.  We  assume  all 
velocity  vectors  lie  in  planes 
through  the  x-axis  and  have  components 
u  and  v  parallel  to  the  x  and  y  axes 
which  are  independent  of  the  angle 
about  the  x-axL3. 


Figure  2.1 

Except  for  an  arbitrary  translation  along  the  y  axis  the  disposition 
of  axes  is  the  same  for  plane  flew. 

With  these  definitions  equations  (2.7)  and  (2.9)  become,  both  for 
plane  and  axial  flow. 


(2.10)- 


|u  ux  +  v  uy  =  -(a"  px)/  Y P, 
I  u  v  ■+  v  v  =  -(a2  pj/  ^  p, 


and 


■where 


(2.11)  /0(upx+vpy)  =  ^p(u/)x+7yfly)5 

u  =  ,  p  =  ,  etc. 

x  dx  *  Fy  ay  * 


The  equation  of  continuity  is 

(2.12)  U/flx  +  T/Oy  +  /°(\.+  Vy+tf  Vy)  “  °> 

where  €  -  0  or  1  for  plane  or  axial  flow. 

Substituting  for  p  ,  p  and  u  p  +  v  /o  from  equation  (2.10)  and 

A  J  A  jT 

(2.12)  into  equation  (2.11)  we  obtain  an  equation 

a2v 

(2.13)  H  ux  4- K(uy  +  vx)  +Lvy  +  <r  “*.0, 


where 


IT  2  2  „  .22 
H  =  a  -  u  ,  K  =  -  u  v,  L  =  a  -  v  , 


independent  of  p  and  /O  ,  This  equation  holds  whether  the  flow  is  ro¬ 
tational  or  not.  In  addition,  if  the  flow  is  irrotational 

(2.14)  Vx  q  =  0, 


or 

(2.15)  v  -  u  =0. 
v  '  x  y 

Shock  Wave 

If  there  is  a  shock-wave  somewhere  in  the  flow,  the  following  equa¬ 
tions  arising  from  the  cpnservation  of  mass,  energy,  rid  momentum  hold 

,.2.2 

across  it  : 


P5  2  n L2  sin2  a  -  (  Y  -1) 

(2.16)  -2= - ~ - W - 


r+  i 


(2.17)  ^2  .  h  8ln 


(Y+  1)  l^2  sin2 

T 


^1  q2  Sin  ^w  “  &)  (Y-  1)  Mj2  sin2  Gw  +  2 
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saXlwtfipau 


I 


Figure  2.2 

In  addition  the  component  of  velocity  parallel  to  the  shock  wave  is  un¬ 
changed. 


(2.18)  qx  cos  ©^  =  q2  ?os  (©w  -  ©). 


In  these  equations  ©  is  the  angle  between  the  velocity  vector  on  the  side 

w 

denoted  tnr  subscript  one  in  front  of  the  shock^wave,  and  the  shock-wave, 
is  q^/n,  and  ©  is  the  angle  between  q^  and  q2  .  From  these  equations 

it  is  possible  to  eliminate  Pg/ph  and  /°  >J (*\  obtaining  the  relation 

“•”>  t  ~  ' 1,1 ' b  - 


An  additional  relation  may  be  read  from  Figure  2.2 

dy  =  ql  ~  *2  ,  d: 

Hx  v_ 


(2.20) 


where  is  the  slope  of  the  shock-wave. 


These  two  equations  are  the  boundary  conditions  which  must  be  satis¬ 
fied  at  the  shock-wave.  Me  shall  call  them  the  shock  conditions. 

Rotation  Introduced  by  the  Shock  Tfave 

From  equations  (2.16)  and  (2.17)  it  is  apparent  that  the  entropy- 
jump  across  the  shock, 


/l^Y 

_pl/ 

w. 

(2.21)  A3  =  s2  -  =  cv  log 

3f  ©  .  Therefore  the  < 

w 

wave  only  if  the  shock  is  straight  in  any  x-y  plane . 


is  a  function  of  © Therefore  the  entropy  is  constant  behind  the  ’•ock- 


Let  the  subscripts  0  and  3  refer  to  stagnation  before  and  behind  the 
shock-wave  if  there  is  one.  Then  from  equations  (2.8)  and  (2.6) 


P//0 


r  _ 


p3//y  > 
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and 


r  P3//»3 


=  1  -  q 


(2.22) 

(2.23) 


P  =  P3(l  “  q2)  '  =  P3  B(q), 

/° r/°3(l  -  q2)5^  =  /®3  A(q), 


where  p0  and  /©_  vary  from  streamline  to  streamline. 

3  '  ^ 

Let  us  attempt  to  define  a  stream  function  z(x,y)  by  the  equations 


(2.24) 


zx  =  -  A  y  v 


and 

(2.25)  zy=Ay’:u, 
where  £  is  as  defined  before. 

Clearly  from  this  definition  V  z  is  perpendicular  to  q.  The  definition 
is  only  justified  however  if 

(2.26)  (zx)y  =  ^ 2y)x  • 

lb  prove  this,  consider  the  continuity  equation 


Differentiating : 


V  .(/O  q  )  =  V  .  (b)  a  q  J  =  0. 


Ad/33  (  V  a.  5)  +  V  •  (A  q)  =  0. 


But  Vz*q  -  0  and  therefore 
V.(A  q)  =  0, 


(A  y  e  u)  +  -iL  (Ay  *v)  =0. 

Referring  tc  equations  (2.24)  and  (2.25)  we  see  that  (2.26)  is  satisfied. 
From  the  equation  (2.2)  we  see  that 


(2.27) 


Pq  /°q 

—  =  ^5-  -  h(«). 


Thus  we  have,  substituting  in  (2.22)  and  (2.23) 

(2.28)  p  =  PQ  h(z)  B(q), 

and 

(2.29)  h(i)*(«). 

These  equations  hold  everywhere  if  we  let  h(z)  =  1  before  the  shook.' 

We  may  now  obtain  a  second  equation  to  replace  equation  (2.15)  when 
the  flow  contains  a  shock-wave,  lb  do  this  we  differentiate  equation 
(2.28)  logarithmically: 


10 


or 


Pv  =  - 


A  (u\  +  rVi 

a 


but  from  equation  (2.10) 

PX=-4(uil  H  U  ), 
p  a  * 

and,  therefore. 


a2  h»  e 

WUVHTT  =  u  u_  +  vu_ 


or  finally 
(2.30) 


„  _  _  £  a  A  h1 

Uy  x  ~  7  ~  1 


This  equation  reduces  to  (2.15)  if  there  is  no  curved  shock-wave  since 
h'(z)  is  then  zero. 

In  order  to  use  equation  (2.30)  it  is  desirable  to  express  -(r-l)h' / 
(rh)  in  terms  of  velocity  components .  This  may  be  done  using  equations 
(2.16),  (2.17),  (2.18),  (2.19),  and  (2.20).  Let  us  write 

^3  =  ^/° 2/  i )  [c/y <°oVc^° 2/ 3) j  /°o  • 

From  equation  (2.29)  1 

£1  *  a  -  ifF*  , 

*0  1 

^  =  (1  -  qj)1^  i 

similarly,  using  equations  (2.17),  (2.6),  (2.19),  and  (2.20) 

/ 0 

2  r+i 

P7‘  7=1 


and 


(tU  -  b) 

7^,’ 


where 


b  =  ^ 


7*1 


JL+  2 


Furthennore,  using  equation  (2.19) 

t  2  .  2  2 

1  -  q2  =  1  -  Uj  -  t2 

=  2qq  (u2 --hc^Uj  -  e) 

-T>TI 
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■where 


e  =  -Tl1  q^  +  ~*r~ 


Therefore 


_  /°3  _  A(q1)(r+1)  r/(r-1) 
'o  Cr-l)(2q1)1^r-1^ 


(u2  -  b) 


Cu2  -  )C  y  u5-e) 


1 1Y 


so  that 


(2.31)  -<•££>  £  «  <“2^) 


dua 


C«2-*)(«2“  u2~e)  dz 

We  shall  call  this  quantity  g(z)  and  rewrite  equation  (2.30) 

(2.32)  vx-„  =4iU 


2.1  J 


We  summarize  these  results  in  equations 


a) 


b) 


2 

H  u  +  K(u  +v  )  +  L  v  +6  =  0 

X  v  y  x'  y  7 

-  2  2  -  t-2  2  _ 

H  =  a  -u  ,  K  =  -uv,  L  =  a  -v  ,  € 


_  f  0  plane  flow 
\l  axial  flow 


V  -u  =£.f  B-g 
x  y  2 


B  =  (1  -  q2)  y'Kr-^) 
^  -  q,  '2 

b«)  g  = 


(u2  -qir 


(uj-bJCuj-  |-  )(  ?"  ^-e) 


dz 


c) 

d) 

e) 


b  s(£t)  k  +  Mr)-V  +(^) 

dz  =  y  6  A(-v  dx  +  u  dy) 

P  =  P0  h  B 

/°  =/°0  h  A 
A  =  (1  -  q2)1/^'-!) 

h  _  A(q1)  (Y* i)rKr-V 

(T^kv17^7^ 


Uj  -b 

T 


L<V  q)(/^)V>] 


>'/(>'- 


mr-1) 
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j  f)  a2=^(l-q2) 

^  g)  v2  (b  -  Ug)  =  (q1  -  Uj)2^  -  d) 

d  =  [(r-i)/(r+i)J  i- 

*  ^ 

h)  a*  =  (qx  -  «2)/v2 

i)  y  =  F(x)  equation  of  contour  of  given  body 

i)  v  =  u  F»  (x) 

k)  z  =  0 

Characteristic  Equations 


Although  it  would  be  possible  to  solve  directly  by  numerical  methods 
equations  2.1a, b,c  subject  to  the  boundary  conditions  2.1g,h,i,j,  and  k 
(being  careful  to  satisfy  the  cause  and  effect  condition  of  Courant  and 
Lewy)  we  have  preferred  to  work  instead  with  characteristic  equations 
to  be  derived  below. 

Hyperbolic  partial  differential  equations  2.4  differ  from  elliptic 
equations  in  that  the  solutions  may  have  derivatives  of  certain  orders 
which  are  discontinuous  across  certain  curves  (if  there  are  two  inde¬ 
pendent  variables)  called  characteristics.  In  the  case  of  supersonic 
flow  these  curves  are  called  Mach  lines  in  honor  of  the  physicist  who 

2  3 

discovered  them  with  shadowgraphs.  It  is  shown  that  if  the  system  is 
of  second  order  there  are  two  characteristics  through  each  point  of  any 
region  where  the  differential  equations  are  hyperbolic;  that  is,  the  flow 

1  7 

supersonic.  Following  the  procedure  of  Frankl  and  Aleksieva  we  shall 
introduce  these  characteristics  as  a  basis  for  a  curvilinear  coordinate 
system. 

Let  oc  (x,y)  =  constant  be  the  equation  of  a  family  of  curves  and 
/3  (x, y)  =  constant  be  the  equation  of  anotner.  Let  OC  and  have  contin¬ 
uous  derivatives  with  respect  to  x  and  y.  If  we  introduce  OC  and  /? 

2  4 

as  independent  variables  ,  assuming  the  Jacobian  oc  and  /Q  with  respect 
to  x  and  y  is  zero  or  infinite  only  at  isolated  points,  equations  2.1a 
and  b  become: 

2.3  R»  Courant  and  D.  Hilbert,  Methoden  der  Mathematischen  Physik,  v.  II. 

2.4  The  method  of  characteristics  as  used  by  3auer  and  flollmien  does  not 
use  coordinates  constant  as  characteristics  but  uses  instead  an  infinity 
of  affine  coordinate  systems  based  on  straight  lines  parallel  to  tangents 
to  two  characteristics. 
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(2.33) 


(2.34) 


(Ho+Kocy)^  +(K«x+  L  Ofy)^  = 


7 

-°Vuo<  +“xT«  =  *V,s  -4xT*  +-Z4JS 


-(H*/»y)V 


If  u  and  v  were  given  continuously  differentiable  functions  of  /$  on 


a  surface  ex  =  constant,  equations  (2.33)  and  (2.34)  could  be  solved  for  u^ 
and  vw  if  and  only  if 


(2.35) 


H  <*  +  K« 

x  y 


Ka  +  Lot 
X 

0 < 


X 


7 


t  0. 


In  this  case  u  and  v^  would  also  be  continuous.  If  u(,0), 
v(^),  Ct(x,y),  /3(x,7)  possess  higher  order  continuous  derivatives, 
it  will  be  seen  by  differentiation  that  the  higher  derivatives  of  u 
and  v,  and  u a  >  v/3  >  etc.  with  respect  to  a  are  determined  by  equa¬ 
tions  (2.33)j  (2.34)?and  their  derivatives.  Therefore,  if  a  =  con¬ 
stant  is  a  characteristic,  equation  (2.35)  must  fail: 

(2.36)  H«x2  +  2K0tx«y  +L«y2  =  0. 

Similarly  it  will  be  found  that  /3  (x,y)  most  also  satisfy  equation  (2.36) 
in  order  to  be  a  characteristic. 


To  be  precise  we  may  define  oc  (x,y)  and  (x,y)  by  the  conditions 


and 


(2.37) 

H«x  +.  (K  -  R)  OC  y  =  0 

(2.38) 

H/3X+  (K  +  R)/3y  =  0 

(2.39) 

J  O  A 

R  =  a  /q  -  a* 

HL 


=  K2  -  R2 


with  some  boundary  conditions  to  be  stated  later. 


From  the  equations 
r 


< 


x<x  ax  +  x/3  =  1  * 

ax +  \  = 0  * 

Xa  ay  +  x/3  =  0  ' 

+  v,  A  =1  j 
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we  find  that 


xa=/3y/A  , 

*6=-VA  * 

(2.39)  7<x  =  ~/V A  » 
y*  -  ax/A 

A  =  ax  /5y  -<*y  /?x  . 

and,  substituting  into  (2.37)  and  (2.38)  we  get 

(2.40)  H  y*  -  (K  +  R)  xa  =  0,  or  (K  -  R)  y^  -  L  x^*  0 
and 

(2.41)  Hy^-(K-R)  X/3  =  0, 

(2.42)  H  U<x  +  (K  -  R)  va  +  xa(£-2-I  +  S^£i)=  o, 

and  2 

(2.43)  (K  -  R)^  +Lv/3  +  J/3  (6$j-  -g^-.g.)  =  0. 

Together  with  equations  (2.24)  and  (2.25)  and  the  boundary  conditions  2.1g, 
h,i,j,  equations  (2.40)-(2.43)  may  be  used  as  a  basis  for  computing  numeri¬ 
cally  plane  or  axial  flows .  If  the  flow  is  plane  and  iirotational,  then 
it  is  preferable  to  introduce  the  velocity  components  as  independent 
variables  because  the  equations  then  become  linear. 

Other  variables  which  are  better  adapted  to  computation  of  certain 
flows  are  q,  the  magnitude  of  the  velocity,  and  0,  the  angle  between  the 

/  2  2  2 

x-axis  and  velocity  vector,  or  /  (q  -  a  )/a  *  p  and  tan  0  =  t;  it  may 
sometimes  be  advantageous,  to  couple  to  z,  a  function  </>  defined  by  the 
equations 

(2.44)  0x=Gu 

(2.45)  <p7  =  Gv 

and  G  must  be  chosen  so  that  )  =  (  ^>  )  .  $  reduces  to  the  velocity 

x  y  7  * 

potential  for  irrotational  flow.  For  future  reference,  we  include  charac¬ 
teristic  equations  in  these  variables  in  our  summary. 
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b)  H  /6X  +  (K+R  )/3y  =  0 

c)  (K-R)  ya  -  L  x^  =  0  or  ya  =  Axa 

L  =  a  -  v 
X=  L/(K-R) 

d)  H  y 6  -(K-R)x^  =  0  or4Jy^  =  x^ 

0)=  h/(K-R) 

e) .  H  ua  +  (K-R)  va  +Xa(  €  tZ  +  ^.11)=  0 

40  V  va  +  xa(£K3?+Eq?)=0 

«  r  a2v  ,,  -  BR  y  g 

P  -  -rr  »  "  — 2 - 


f) 

(K-R) 

>  +  L  > 

+  ^  (‘TT- 

BR_ 

yggy_ 

V+At4  +y/3 

(*  IPS  ‘ 

Q 

K=K 

)  =o 

g) 

a  y£ 

- 

/  2  2 
G  Z/3  /q  -  a 

=  0 

g') 

*V 

Ou,  <p7  ■■ 

=  Gv 

h) 

a  y  ^ 

A  +“ 

/  2  2 
z«  /  q  - a  : 

=  0 

h«) 

zx  = 

-  yf  Av, 

ZY  ~  jC  Al 

a 

i)  q  /sLl^  +  eflf+-  r€5i^+L_^L  ; 

a)  -o  i£ST  + Csis£-ia!/i 


i)  q  yH.ii  +  e 

ct  aq  < 


a  cos  9 


/  2  2 

)  n)  y a  (t  +  p)  =  (tp  -  1),  t  =  tan  Q,  p 


a 


°)  (p  -  t)  =  (tp  +  i) 

#2  2 

p)  dz  =  —  (dy  -  tdx),  q2  =  S  ■  £  2 

'  JZ  1  +  q  p 


1+t 


*2  _  r  -l 
q  - 


r+i 


2  *2 

q)  (1  +  t2)  f(p)  Pa  +  ta  +  |  (y^-  txa)  =  0,  P2cV“  q 


y  wet 

r)  -(1  +  t2)  f(p)  Pfl  +  +  |  (3jg  -tx^  )  =  0 

f(P)  =-il--.3l.P . 

o  %  .'<2  2 

(1+P2)(1-K1  P^) 


r (1  -  q  ) 


«,  r  <P_cav 

S  '  TPE  ~  "'y^-g?  ’ 


ip  r  Q  _  BR  y*  g 
”  2(K-rj 


Section  3»  Numerical  Solution  of  Boundary  Value  Problem 
Introduction 


Typical  of  the  characteristic  equations  which  may  be  used  to  com¬ 
pute  the  supersonic  steady  frictionle3s  flow  past  a  given  olane  or  axi¬ 
symmetric  body  are  equations  2.II  c,d,e,f,  and  2.1c  with  boundary  con¬ 
ditions  2.1  g,  h,i,j.  Accordingly  we  shall  describe  the  procedure  we 
use  in  terms  of  these  equations,  This  is  no  restriction,  of  course,  and 
any  of  the  other  Sets  of  equations  may  be  used,  e.g.,  2.IIn,  o,  p,  q,  r 
with  boundary  condition  deduced  from  2.1g,h,i,j» 

We  shall  consider  the  ca33  of  axisymmetric  rotational  flow;  the 
ca3e3  of  axisymmetric  potential  flow,  plane  rotational  or  plane  potential 
flow  may  be  treated  the  3ame  way  with  several  simplifications. 

Consider  then  a  supersonic  flow,  uniform  at  infinity,  past  a  body 
of  revolution  ABDEF  (Figure  3.1).  The  fact  that  we  have  assumed  the 
flow  supersonic  implies  a  restriction  relating  the  free  stream  Mach 
number  and  nose  angle  of  the  body.  Indeed  if  the  nose  is  blunt,  or  if 
it  i3  pointed  but  the  free  stream  Mach  number  i3  less  than  some  number 
greater  than  one,  it  is  known  that  the  3hock  wave  is  detached  from  the 
body  and  crosses  the  axis  normal  to  it.  from  equation  (2.17)  it  fol¬ 
lows  that  q»  would.be  less  than  [(Y  -l)/(  T'+l)!  i.e.  subsonic. 

For  cones  there  i3  a  half-vertex  angle  (about  52°  34'  in  air)  above 
which  the  3hock-wave  is  detached  at  all  Mach  numbers.  For  each  smaller 
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cone-angle  there  is  a  Mach  number  above  which  the  shock  is  attached  to 
the  vertex  and  there  are  two  conical  shock-waves -each  of  which  corresponds 
to  a  mathematical  solution  of  the  flow  problem. 

Intuitively  it  is  clear  that  for  other  bodies  of  revolution  than 
cones  there  exist  one  or  more  solutions  (probably  two)  with  attached  shock 
if  the  nose  angle  is  small  enough,  the  curvature  negative  or  zero,  and 
the  Mach  number  large  enough.  As  far  as  we  know,  this  has  not  been  demon¬ 
strated  mathematically,  although  the  paper  by  Frank!  and  Aleksieva1, '  con¬ 
tains  a  theorem  which  the  authors  believe  could  be  extended  to  do  it.  We 
shall  assume  that  it  is  true  and  that  we  have  this  case  before  us.  We 
shall  assume  furthermore  that  the  solution  is  a  continuous  function  of  the 
boundary  y  =  F(x)  in  the  sense  that  if  we  replace  a  small  section  of  the 
nose  by  a  straight  line  A'B  tangent  to  it  at  the  point  B  of  juncture  and 


Figure  3«1  Figure  3*2 


let  P  approach  A,  the  flow  about  A'BDEF  will  approach  the  flow  about 
ABDEP  i «e . 

lim  x'  (  a ,  /3  )  =  x(  a  ,  /3  )  , 

B  — *>A 

lim  u'  ( a,/3 )  =u(a,/6) 

B— A 


etc. 

In  summary  we  assume  that  the  given  body  of  revolution  has  a  con¬ 
tour  characterized  by  the  equation  y  =  F(x);  that  F*  exists  except  at 
isolated  points,  and  is  less  than  or  equal  to  zero  everywhere  between 
x  =  0  and  x  =  Xj,  ,  except  at  a  finite  number  of  points  such  as  D  and  E 

where  F'(x)  may  be  discontinuous;  that  F' (0)  «tan  52°;  and  that  is 

large  enough  so  that  the  shock-wave  is  attached. 

3.1  (j.  T.  'Taylor  and  J.  W.  Maccoll,  Proc.  Roy.  Soc.  of  London,  Series  A, 
vol  139,  1933;  p.  278-299. 
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Approximations  at  the  Nose 


At  the  nose,  in  accordance  with  an  assumption  stated  above  we  re¬ 
place  the  contour  to  the  left  of  B  by  a  straight  line  tangent  to  the 
given  contour  at  B. 

3  2 

Since  the  flow  at  a  point  P  is  independent^*  of  changes  made  in 
the  region  bounded  by  the  two  half  characteristics  farthest  from  the 
velocity  vector  at  P,  the  flow  in  the  region  A'BC,  bounded  by  A'B.  a 
characteristic  BG,  and  the  shock  A'C,  is  precisely  Taylor-Maccoll3*l 
flow  over  a  cone.  That  is,  u  and  v  are  constant  on  lines 


(3.1) 


x  -  x 


=  t 


“A1 


through  the  nose  A’.  We  may  therefore  seek  immediately  the  values  of 
x,t,u,v,z  for  equally  spaced  values  of  some  variable  such  as  y  along 
the  characteristic  BC.  The  differential  equations  for  u,v,and  t  may 
be  found  more  readily  from  equation  2.1a,b,  2.IId  than  from  Taylor  and 
Mac coil's  equations.  Since  u  and  v  depend  snly  on  t, 


_  du  ,  _  du/dy  t  _  u' 

ux  “  M  \  -  y - V 


cRT  x  dt/dy  y 
where  we  denote  by  a  prime,  d/dy.  Similarly 


u  = 

y 


Vx 


u' 

v' 

V 


and 


v  = 
v 


V' 

V 


£ 

y 

t 

y 


£ 

y 


Therefore  from  2.1a, b  (g  is  zero  in  A'BC  since  the  shock  is  straight) 


u'(K  -  tH)  +  v'(L  -  tK)  + 


2  ,  . 
a  vt' 


=  0, 


u'  +  v't  =  0. 

In  addition  it  follows  from  equations  (3*1)  and  2.IId  that 


(3.2) 


* .  _  t(K  -  R  -  tH) 
t - y(K"-"R7  ■ 


Using  equation  (3*2),  the  equations  above  (3.2)  may  be  solved  for  u'  and 
v' : 


(3.3) 


u'  = 


2  + 
-a  vt 


y  £t(K  _  fc)  -Lj 


5.2  R.  Courant  and  D.  Hilbert,  Methoden  der  Mathematischen  Physik, 

voi.  n,  p.  307  . . . 
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(3*4)  v' =  fm*  -  SR}  * 

To  these  equations  vre  may  add  three  initial  conditions 

(3.5)  t  _  yB 

0  *B-^. _  ' 

(3.6)  uo  =  qo/  /l  +  to2  , 

and 

(3.-7)  Vq  =  qoto//l  +  to2  . 

Integration  is  to  stop  when  u,  v  and  t  satisfy  equation 

(3.8)  (y+  l)t  (u+vt)(tu-v)  =  (/-l)  £l  +  t2  -  (u  +  vt)2J 

derived  from  shock  conditions  2.1  g  and  h  by  elimination  of  q.. .  q.. 
is  then  determined  from  equation 

(3.9)  qx  =  u  +  vt. 

The  value  of  q^  so  obtained  will  vary  with  qQ.  Therefore  in  order  to 
obtain  the  flow  at  a  prespecified  Mach  number  or  value  of  q^  we  shall  have 
to  modify  qQ  in  a  way  governed  by  the  variation  of  q^» 


Distribvfting  the  data  at  equal  intervals  of  y  yields  a  poor  distri¬ 
bution  in  the  hodograph  plane.  It  is  better  t^  use  something  dependent 
on  the  velocity  as  independent  variable,  e.g.  -  *  T.  If  this  be  done, 
it  may  be  shown  as  above  that  the  differential11  equations  are 


r 


u» 

II 

gfs 

II 

-yu 

x+yT 

z' 

=  -  y  a 

A  f(va  +  Uy/q2  ■ 

X' 

t  2 

*  (a  - 

u2)f 

y 

=  -  (uv 

+  a\/q2  -  a2)  f 

V 

=  uT, 

f= 


a  v(x+yr; 


with  initial  conditions 

(3.5')  10  =yB/CIB-IA,)  , 

U  =  q  /  J 1  +  I  3  , 


z  =  0  , 
o  * 


and 


Xo  =XB  * 
*Vo  =  yB  J 
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and  terminal  conditions  to  detemine  q  and  q.  : 

o  1 

(3.6')  [(y2+x2)u  ~  x(ux^vy)J  =  (y^)x[x2+y2  -  (ux+vy)2J 

_  _  ux  +  vy 

ql  "  x -  * 

If  the  problem  is  to  be  done  by  hand,  equations  (3.2) ,  (3.3),  and 
(3*4)  may  be  solved  first;  r.  and  /I  found  later.  With  the  ENIAC,  how¬ 
ever,  this  would  waste  time  since  more  cards  would  have  to  be  read. 

From  equation  2.1c  we  find  that 

(3.10)  z»  +  [u(K-R)  -  vh]  =  2|£  [a  u  -  v /q2  -  a2]. 


Freedom  in  Choice  of  OC  and  /3  . 

As  for  /3  ,  let  us  consider  more  generally  the  detennination  of  oc 
and  /3  throughout  the  whole  region  in  which  we  shall  seek  the  values  of 
x,  y,  u,  v,  z.  Suppose  that  we  assign  /3  arbitrarily  on  BC  increasing 
from  B  to  C  (Figure  3*3)*  Assuming,  as  we  have,  that  the  characteristics 
have  no  envelopes  and  that  therefore  there  is  one  characteristic  of  each 
family  through  each  point,  the  values  of  /3  in  BC  determine  the  values 
of  through  the  region  BCD  but  have  no  effect  on  the  values  of  a  in 
that  region.  Accordingly  we  are  free  to  make  oc  an  arbitrary  increas¬ 
ing  function  on  BD.  This  will 
detemine  ot  in  the  region  BBEC 
but  will  not  affect  the  values 
of  in  CDE.  (E  may  be  at 
infinity  since  BE  may  not 
intersect  the  shock  wave.) 

Step  by  step  it  is  seen 
that  /3  may  be  assigned 
arbitrarily  (we  shall  make 
it  increase  from  B  to  H) 
along  BCEH  and  CX  my  be 
assigned  arbitrarily  along 
BDFI  (we  shall  make  it 
increase  from  D  to  I). 

Indeed  we  may  make  oc  = 

0i  ( fl)> 

~  61 {oc)  along  EDFI 

and  -fa  (oc),  « 
along  CEH  where 


Figtcre  3*3 

d?  (fi)  and  (ot  )  are  non-decreasing  functions. 

In  order  to  make  the  map  of  the  region  BDFIHECB  on  to  a  portion  of  the 
OC/S  plane  a  one  to  one  map,  it  is  necessary  tc  forbid  the  maps  of  3DFI 
and  CEH  to  have  any  points  in  common.  Because  there  are  two  arbitrary 
functions  at  ora's  disposal  in  the  assignment  of  parameters  there  are  many 
choices  available .  One  could  .let  -  y  on  BCEH  and  oc  -  x  on  BDFI ;  one 
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could  let  oi  =  /3  =  x  on  BDFI  or  on  GEHG  (in  the  latter  case,  this  leaves 
cx.  undefined  downstream  of  the  characteristic  oc  =  o«  which  may  intersect 
the  body  somewhere  between  D  and  I);  one  could  let  -  y  on  CEKG,  or  if 
s  is  arc  length  along  EDFI,  s  >  is  arc  length  along  CEHG  and  O'  is  curva¬ 
ture  along  either,  one  could  let 

^  =«q  +«x  <r(s)  on  EDFI  , 

=<3  +A  o^s*)  on  CEHG 
s'  '  o  'i 

where  ocQ, /3n,  are  positive  constants,  so  computations  carried  out 

on  a  square  grid  in  the  Oi  ,  plane  would  correspond  in  the  physical 
plane  to  a  fine  grid  where  the  velocity  is  changing  rapidly  and  a  large 
grid  where  it  is  changing  slowly.  Dther  choices  will  suggest  themselves 
to  the  reader.  To  facilitate  the  simultaneous  programming  of  whole 
families  of  flews  for  ENIiH  and  other  machines  we  shall  in  this  report 
usually  map  the  shock-wave  onto  a  line  of  slope  one  and  the  given  con- 
tour  onto  lines  of  slope  \  or  1. 

Corners 

From  a  point  on  the  given  contour  such  as  D,  F  or  I  on  Figure  3*4 
where  the  slope  is  decreasing  but  discontinuous,  it  is  known  that  a  family 
of  characteristics  OC  =  constant  emanates,  but  only  one  characteristic 
=  constant.  Therefore  points  D,  F  and  I  must  map  onto  horizontal 
line  segments  D'...D1,  F1 ....F1  and  I»...I«.  Once  again  the  assignment 
of  oc  is  arbitrary. 


Terminal  Boundary  to  Flow  Computations 


Since  the  velocity  at  G  is  independent  of  the  velocity  at  points  to 
the  right  of  characteristic  HG,  there  is  no  reason  to  compute  the  flow 
beyond  this  curve  if,  for  example,  the  pressure  distribution  along  A'BDFIG 
is  desired.  The  flow  can  be  computed  in  the  region  HGF  bounded  by  the 
characteristic  HG.  It  will  have  to  be  computed  also  in  some  region  HGLK 
in  order  to  determine  the  flow  in  the  base  region.  At  the  present,  there 
is  no  known  satisfactory  way  of  doing  this  because  the  non-viscous  steady 
flow  model  is  not  applicable  to  the  wake,  as  a  glance  at  a  shadowgraph 
of  a  wake  will  suggest.  Thus  the  object  of  this  report  will  be  to  deter¬ 
mine  the  velocity  at  a  net-work  of  points  in  the  region  R  or  R'  shown  in 
Figures  3»4  and  3.5. 


Returning  to  the  determination  of  the  flow  along  the  characteristic 
BC,  we  shall  simply  set  ot  =  O  and  /?  =  /3  (y  -  y  ).  Summarizing:  along  BC 

O  D 


a) 

b) 

c) 

d) 


dt  _  t(K-R-tH) 
3y  y  (Y  -  t) 


=  _£. 


x-x 


du 

Hy 

dv 


A' 

-a2  v  t 


y  [t(K-R)-Lj 
a2v 


y  [  t(K-R)-L  J 

e)  a§  -  &  0<K-R)-*h] 


f) 

g) 

h) 


#  =  /?0Cy  -  yB) 


oc  = 


7i 


t  -  B 


u  = 


A' 


v  = 
0 


Vc 


1+t 


A  + 12 

0 

D  O' +«V'W.)f Vv^.)  “O'-1)  [Kh«,«,V!] 

x  +  yT) 

,  A  2T 

+  u  yq  -  a  ) 


i) 

(/ 

+1)VWi 

)) 

ql 

=  \  +  V. 

k) 

!X> 

*  ar  “  -  ' 

1) 

z* 

=  -ya  Af(v  a 

m) 

x» 

"  (a2  -  u2)f 

n) 

r* 

=  -  (uv  +  a 
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jo)  v  =  u  T 

p)  f  =  -  y  u 


3.IS 


|{l/(K-R)}  X  -  yj/[a2v(x  +  yT)J 


q)  q2  =  u2  +  v2,  a2  =  — (1  -  q2) 


2  2 


2  2 


r)  H  =  a  -u  ,  K  =  -uv,  L  =  a  -v  ,  R 

s)  A  =  (1  -  q2) 


/  2  2 
=  a  /  q  -a 


To  =  V(XB  “  XA'5 

»  =  <l/ 


t)  <  z0  =  0 


Xo  =XB 


yo  =  yB 


u)  /5  =  (T  -  T0 J/CTLj^  -  To) 

[2  2 
(y  +x  )u  -  x(ux+vy) 

w)  qx  =  (ux  +  vy)/x 


)1  -/V"A  x  r 

J  \T+T j  [ 


2  2 

X  +y  -  (ux+vy) 


2] 


t,  u,  v,  2,  x,  /d  are  determined  at  equally  spaced  intervals  of  y  by 
solving  equations  3.1  a,b,c,d,e,f  with  initial  conditions  3.1h  and  termi¬ 
nal  conditions  3*1  ij jj  or  u,z,x,y,v,  /3  are  determined  at  equally  spaced 
intervals  of  T  by  solving  equations  3*1  k — n,  u,  with  initial  conditions 
3. It  and  terminal  conditions,  3*1  v,w. 


Contour  Process 


Once  the  initial  data  has  been  determined,  flow  variables  may  be 
found  at  the  intersection  of  the  given  contour  and  a  characteristic 

,/3  =  constant  through  a  point  on  BC.  The  process  involved  in  doing 
this  we 


shall  call  the  contour  process.  It  will  be  used  every  time  the  flow 
variables  are  to  be  found  at  a  new  point  on  the  contour. 
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Since  only  oc varies  in  passing  from  B,  to  B..  two  of  the  equations 
to  be  used  are  2. lie  and  2. lie 

(3»U)  y  oc  ~  ^  ^  oc  * 

and 

(3.12)  &>ua  +  Vot  +  ^C£  P  +  Q)  =  0. 

Actually  Q  nail  be  zero  since  the  shock  is  straight  near  z  *  0.  The 
other  three  equations  are  the  boundary  conditions  2.1  i,j 

(3.13)  y  =  F(x), 

(3.14)  v  =  uF'(x)  =  uG(x) 

and  the  condition 

(3.15)  z  =  0. 


A  simple  procedure  which  can  be  used  to  find  the  flow  variables  at 


is  simply  to  replace  by 


ai^be_^se' 

A  oc 


etc.,  obtaining  linear  equa¬ 


tions  for  x,  y,  u,  v.  The  justification  for  this  procedure  is  given  in 
Frank!  and  Aleksieva's  paper  which  applies  to  our  problem  once  we 
assume  that  the  flow  on  BC  is  correctly  determined.  Instead  of  using 
equation  (3-13)  as  it  is  we  prefer  to  differentiate  it  in  order  to 
obtain  a'  linear  equation  in  xu  and  y_ 

"l  U1 

(3.13')  g  =  G(x). 

Let  us  denote  the  map  of  B^,  B  and  in  the  &-  plane  by  /  ,  a, 
and  no  name,  (Figure  3*7)  and  denote  the  corresponding  x,y,u,v,z  by  xw  ,y#  , 
U/  *  Y  *  *  V  V  V  V  V  X>  7i  u*  v>  z’  111611  we  ^t,  using 

the  suggested  procedure,  the  following  linear  equations 


(3.16) 

y-ya  =  F'  (xa)(x  - 

(3.17) 

y-y Xje  (x-x/)' 

(3.18) 

v  =  u  G  , 

(3.19) 

(u  -  Uj)  +  V  -  V^  +£■■ 

(3.20) 

z  =  0  j 

.,  solved  for  x,  y,  u,  v,  z  give 

(3.21) 

3C  =  [yz  -  ya  +  F'4  xa  -  X 

(3.22) 

7  =  Yg  F^  -^ya  + 

(3.23) 

N 

II 

o 

j e 

^7 


(X-X£)  =  o. 
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/ 


and 


(3.24) 

(3.25) 

v  =  G  fay 

T 


4jj- (X  -  Xj  ) J  /(CJjf  +  G), 


j-4e -  (x  '  *£  l]  /(“/+  °)* 

S  S 

This  represents  the  simplest  possible  contour  process.  Clearly  many 
refinements  may  be  made.  For  example,  it  would  be  preferable  to  evaluate 
X  ,  CO  and  P  at  m,  the  midpoint  of  A  and  the  desired  point,  and  G  at  M, 
the  midpoint  of  a  and  the  desired  point.  This  may  be  done  by  extrapola¬ 
tion  or  by  integration.  We  shall  reserve  discussion  of  such  refinements 
for  the  next  section. 


General  Process 

Once  the  flow  variables  have  been  found  at  or  any  other  point  on 

the  contour,  the  next  step  is  to  find  them  at  a  point  P  at  the  interetto- 
tion  of  a  characteristic  p  ~  constant  through  B2  on  the  initial  line  and 

a  characteristic  oi  =  constant  through  D^.  ISore  generally,  given  the  flow 

variables  at  any  two  points  j^and  u  not  on  the  3ame  characteristic  we  may 
find  them  at  the  intersection  P  of  characteristics/^  =  constant  and  ot  = 
constant  through  JL  and  u  respectively-  We  3 hall  call  the  process  for 


doing  this  the  general  process. 

The  equations  to  be  used  to  this  end  are  equations  2. II  cdefh* 

(3.26)  =  ^  x  *  » 

(3.27)  aiy^  =  Xg  , 

(3.28)  oi u^  +  va  +  xa  +  A)  =  0  » 

(3.29)  Uj+  Xvg  +  y^  (e-^R  -  £!£)=  °, 

and 

(3*30)  dz  =  yA  (-  v  dx  +  u  dy). 

The  method  again  is  to  replace  partial  derivatives  by  difference  quotients 
ttiare  is  no  need  for  the  grid  sizes  in  a  and  ft  to  be  equal  since  the 
difference  equations  do  not  contain  A  Ot-  or  A/3  : 
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(3.31)  y  -  X  x  =  y^  -  X* 
(3-32)  wy-x=a>yu-xu  , 


(3*33)  to  u  +  v  =  tou^  +  -  (x  -  )(S  +  f)  , 

(3.34)  u  +  \v  =  uu  +  ”Xvu  -  (y  -  yu)  (S  -  T)  , 

and 

(3.35)  2z  =  +  zu  +  y  A  [-(2x  -  x^  -  xu)v  + 

where  s=FlandT=rS* 

Instead  of  evaluating  the  quarrtitites  a  ,  y,  I,  etc.  at  the  midpoints 
of  JL  and  P,  and  u  and  P  we  suggest  in  this  simple  general  process  evaluating 
them  at  the  midpoint  of  JL  and  uj  i.e.. 


(2y-  jt  -yu)  5]  , 


I  =  Afq2), 

etc.  This  simplifies  the  computation  and  is  all  that  is  justified  until 
a  more  careful  procedure  is  describe^, in  the  next  section.  Solution  of 
equations  3.31-3.35  yields  x,  y,  u,  v, z  at  P. 

(3.36)  x  a  |j(y^  -  Ax^  )  -  (i>yu  -  xu)J  /(I  -%u>  ), 

(3.37)  y  =  j(y£  -  Xx^  )  -  X  (wyu  -  xu)J  /(l  -  Xw), 

(3.38)  u  =  [{uu+  X  vu-  (y-yu)(§-T)}  -  X{^  -(x-3^)(S+T)]1  /  ' 

(1-  Ju>), 

(3.39)  v  =  |{«^+  ^  )  (s+fjj  -  w-[uu+  A  vu  -  (y-yu)(s-T)]I  / 

(1-  Aw)> 

(3.40)  z  =  +  zu  +  yJ  {-(2X-0CJ  -X^)v  +  (2y  -  y^  -  yu)  ujJ  . 

Any  hand  computation  must  be  accompanied  by  checks .  The  values  of 
H,  K,  L  and  R  may  be  checked  by  the  identity 

HL  =  K2  -  R2  J 

the  solutions  (3*36)  -  (3*40)  by  substitution  into  some  of  equations 
(3.31)-(3*34)  the  value  of  z  by  the  formula 

z  =  z£  +  y  I  (x  -  *£  )  v  +  (y  -  y,  )u] 

etc. 
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Shock  Process 


Similar  to  the  contour  process  though  more  complicated  is  the  other 
boundary  process  which  gives  the  flow  variables  at  a  new  point  on  the  shock 
wave.  In  this  case  the  flow  variables  are  known  at  a  point  a  oh  the  3hock 
wave  and  another  point  u  lying  on  a  characteristic /d  =  constant  through  a. 
The  point  P  lies  on  the  shock  wave  and  a  characteristic  (X  =  constant  through 
u. 


Figure  3.10  Figure  3.11 


The  equations  we  shall  use  are  2.II  d,f,  h*  and  2.1  g,h,b* 


(3.41) 

(3-42) 

(3.43) 

(3.44) 

(3.45) 
(3-46) 


Hx  v  * 

u/3  +  A  V/3  +  y/3  {  =  °» 

v2(b  -  u)  =  (qx  -  u)2(u  -  d),  or  ^  =  v«(u,v). 


N 

dz  =  yA  (-  v  dx  +  u  dy). 


The  procedure  is  as  before;  the  coefficients  in  equations  (3*41)  and 
(3.43)  may  be  evaluated  at  u,  the  other  at  a.  The  results  for  x,  y,  u, 
v,  z  are 


(3.47)  =c  =  [ouhaxa  -  xu  +  «u( VJTa )]  /(“A  -  1). 

(3.48)  y  =  [haVu-  ya-ha(\  "  v]  /(»„•>,  “  1)' 

(3.49)  u  =[Aufaua  +  uu  +  Au(vu  -  Ta)  -  (y  -  ru)<Su-V] 

(3.50)  r  =  [faVu  +  7a  +  fa(uu  ‘  V  ‘  Mifa  +  1'>’ 

(3.51)  .  =  *„  +  ya  Aa  [-va(x-xa)+ua  (y  -  ya)]  , 


ha  *  (ql  -  V^a 

fa  '  [^l  -  S)  {  <V%>  +  2<dM1a>}  wa2]  /{X*1  +  Va>]  ' 


where 


When  the  flow  variables  are  determined  at  P,  g  may  be  computed  at 

the  midpoint  of  aP  by  equation  (3.46)  using 

,  u 

du  _  u  -  a 

35s  z  -  z  * 

a 


Corner  Process 

We  have  seen  that  a  corner  such  as  D  (Figure  3.2)  maps  into  a  line- 
segment  D — D  in  the  oc  plane.  Along  DD  u,  v,  and  ex.  vary  although  x, 
y,z  are  fixed.  Oi  may  be  related  to  u  and  v  in  any  practical  way  such  as 


(3.52)  <X  = 


&1-*2>\V2  v  ,  “2V2  '  aiT2aX 

V1  "2  '  v2  “1  5  V2  -  V2 


a  linear  function  of  the  tangent  of  the  velocity  inclination. 


Figure  3.12  Figure  3.13 


Since  x  is  constant  at  D, 

XOi  ~  0  > 

and  thus  from  2. He 

(3.53)  *>ua-  +  v*  =  0  , 

or 

(3*54).  wdu  +  dv  »  0  . 

Ibis  is  an  ordinary  differential  equation  in  u  and  v  whose  solution  is 
an  epicycloid  in  the  u,v  plane.  The  solution  has  been  tabulated  in 
3  3 

various  places^  .  Thus  it  is  simple  to  obtain  a  set  of  points  along 
D — D.  However,  if  more  accuracy  is  required  than  three  significant 
figures,  and  if  it  is  desirable  to  space  the  data  at  equal  increments 
in  v/u,  then  it  is  simpler  to  modify  equation  (3*53)  introducing  v/u  = 
t  as  independent  variable 


S*g»,  N.A.C.A.  Technical  Note  1428,  Dec.  1947, 
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29 


<3-55)  £=imr  • 

and  solve  directly. 


Section  4.  Refined  Numerical  Methods 
Introduction 

In  section  3  the  problem  of  determining  numerically  the  supersonic 
steady  frictionless  flow  past  a  body  of  revolution  was  solved  in  a  simple 
way.  Having  assumed  that  the  flow  near  the  nose  could  be  reasonably  well 
approximated  by  the  Ihylor-Maccoll  flow  with  the  smaller  shock  wave  angle, 

1  7 

the  paper  of  Frankl  and  Aleksieva  *  proves  that  there  exists  a  unique 
solution  and  that  the  approximation  outlined  in  the  last  section  con¬ 
verges  to  this  solution  as  the  mesh  size  approaches  zero  (if  we  ignore 
round-off  errors). 

However,  in  the  last  section,  little  attention  was  paid  to  the  prob¬ 
lem  of  getting  the  best  approximation  for  a  given  amount  of  labor  by  hand 
or  by  machine,  ttiis  is  a  problem  which  can  never  be  completely  solved. 
Nevertheless,  in  this  section  we  shall  examine  some  aspects  of  the  prob¬ 
lem  and  show  how  some  of  the  computations  previously  described  can  be 
done  to  a  given  accuracy  more  easily.  Tfe  shall  indicate  in  several 
places  how  a  problem  should  be  treated  depending  on  whether  it  is  done 
by  hand,  or  by  a  fhst  machine  of  small  memory,  or  by  a  fast  machine  of 
large  memory. 

Systems  of  Ordinary  Differential  Equations  with  Initial  Conditions 


Given  a  system 

(4.1)  7±  -  x),  i  =  1,  2,  - ,m,  y±  =  y±(x), 

of  ordinary  differential  equations  with  initial  conditions 

(4.2)  y±(*0)  ~  y±  • 

o 

The  numerical  method  of  solution  most  commonly  used  for  hand  compu¬ 
tation  at  the  BRL  is  due  to  1'oulton^*^,  although  it  differs  only  a  little 
from  a  method  used  earlier  by  Adams.  Like  most  numerical  methods  it 
assumes  that  the  solutions  may  be  closely  approximated  on  short  intervals 
by  polynomials  of  suitable  degree.  The  polynomial  of  degree  n  passing 
through  n+1  points  has  been  found  in  various  forms  by  Gregory  (1670), 
Newton  (1687),  Waring  (1779),  Lagrange  (1795).  If  Gregory's  formula  for 

(4.3)  fi(3i,....yn>  x)  =  Fi  (x) 

4.1  Bennett,  lElne,  and  Bateman,  "Numerical  Integration  of  Differential 

Equations";  Bulletin  of  the  National  Research  Council  No.  92,  Nov. 
1933,  pp.  747757^ 
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be  integrated  from  to  xn+1  the  following  formulae  result: 

J^nfl  _  -\  o 

Fi(x)dx  =  h  Ui>nfl  -iV  Fi<n+1  -  IjV2  Fi>n+1  ] 

n  u 

+  h  [-  h:  V  3  V*  -  *  4  Fi,n+1  -mV5  FijWl] 

+  h  [-  EOTJ  V6  Fiin+1  -  ••••3  * 


where 


(4.5) 


h  =  x  ,,  -  x 
n+1  n 


Fi,n+1  *  V^ml*  y2,n+l»  ***  ^n+l*  xnM^ 


^ Fi,n-fl  “  Fi,n+1  “  Fi,n 


A,n«  *  V  [VFi,n+l  -  *Fi,n.>  Fi,n+l-2Fi,n+Fi,„-l 


ly3  Fi,wl  *  V[V  3-1  Fi,n+X  -  VJ_1  F1>n3  • 

Moulton's  method  ir  the  following:  values  of  y^  Fi  y ^  2, 

F.  y.  F.  .  are  obtained  corresponding  to  xi  x0  x.  by 

special  means  such  as  Taylor's  series.  These  aie  arranged  on  a  computing 
form  as  follows: 
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A  trial  value  for  F.  .  _  is  then  secured  by  extrapolation  or  otherwise, 

i,  j+i 


and  y.  .  ,  is  computed  from  the  formula 
1,  j+-i- 


(4.7) 


p+ 1 

h,i+i=y±,i  +  \  F^> 

•S  v 


using  equation  (4.4). 

A  corrected  value  of  F 


i>j+l 


is  then  found  from  equation  (4.1)  and 


the  process  repeated  until  no  change  occurs.  The  values  of  x.  7-  . 

i-1  i-L.  1»J+± 

F.  .....  F,  v?  r  -.-i  are  entered  in  the  table.  IfiLues 

ljJ+l*  v  i>J+l  v  m,3+l 

of  J-i  -jo*  7.  -.u  etc.,  are  computed  in  the  same  way  until  the  eni  of 

the  interval  is  reached.  Practically,  differences  beyond  the  third  are 
not  often  used  so  that  the  solution  is  approximated  by  a  fourth  degree 
polynomial.  We  shall  refer  to  it  then  as  a  fourth  order  method.  More 
generally  we  shall  say  that  if  j-1  differences  are  used  the  method  is 
a  jth  order  method,  because  if  the  functions  y^were  anal^ic  we  should 

be  making  truncation  errors  of  the  form  [d^^y/dx^J  ^^+^/(j+l)l  J 
in  each  interval.  It  is  customary  to  choose  j  conveniently  and  then  to 
adjust  the  size  of  the  interval  so  that  xj  F.  .  are  negligible.  This 

may  require  several  changes  of  interval  size  during  the  computation. 

It  will  be  noted  that  if  it  be  necessary  to  halve  the  interval,  then 
auxilidiy  points  must  be  inserted  using  interpolation  formulae.  As  a 
guard  against  errors,  hand  computations  are  frequently  checked  using 
equation  (4.7)  over  larger  intervals  and,  for  example,  Simpson's  Rule 
and/or  Weddle's  Formula. 

Moulton's  method,  though  convenient  for  hand  computation,  has  three 
defects  for  machine  computation  which  are  avoided  by  other  methods.  These 
are:  (a)  the  number  of  quantities  which  must  be  remembered  in  going  from 
x^  to  x^+1,  namely,  |(j+l)m+lj  j  (b)  the  fact  that  early  steps  are 

different  from  later  steps;  and  (c)  the  necessity  of  using  interpolation 
formulae  in  reducing  the  interval  size.  In  setting  up  a  problem  for  a 
machine,  if  there  are  N  registers  available  to  remember  numbers,  then 
l+(j+l)m  must  be  less  than  or  equal  to  N.  Ihus  if  a  fourth  order  approxi¬ 
mation  is  to  be  used,  then  m  must  be  less  than  or  equal  to  (N-l)/5.  For 
example,  N  will  be  about  twelve  to  fourteen  for  the  ENIAC  if  ten  figures 
are  to  be  used.  Let  us  in  future  discussions  let  N  =  13;  thus  only  a 
second  order  system  could  be  handled  by  Moulton’s  method  on  the  ENIAC. 

Most  of  the  other  common  methods  such  as  those  of  Adams,  Steffensen 
and  Milne  are  subject  to  the  same  objections  for  machine  work  and  accord¬ 
ingly  will  not  often  be  used  for  high-speed  machines. 

The  best  known  methods  which  are  less  open  to  objections  (a),  (b), 
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/  o 

and  (c)  are  those  of  Runge  and  Kutta*'  .  The  step  from  x^  to  is 

the  same  for  all  k  with  these  methods .  For  example,  their  fourth  order 
method  is  described  almost  completely  by  the  formulae 

7i,nfl  =  7i,n  +  ^ki,l  +  2ki,2  +  21i,3  +  ki,4^6’ 
where 

ki,l  =  h  fi  'V  »i,n’  —  \n> 
k.;2  =  h  t±  (xn  +  ih,  y1<n  +  ik^,  ....  7m<n  +  ikB<1) 

\,3  =  hfit\*  ih'  Jl,  n  +  V,n  + 

ki,4  =  h  fi  (xn  +  h-  »l,n  +  h,3‘  — >  +  km,3) 

V 

(i  =  1,  2,  ••••,  m) » 

This  method  is  free  from  defects  (b)  and  (c).  As  for  (a)  there  is  one 
stage  of  the  computation  when  xq,  7l,n’  ***’  7m,n'  V  7l,n  +  (ki,l+2kL,2)/6’ 

72,n  +  ^,1  +  2k2,2^6’  7m,n  +  <\,1  +  2km,2^6,  ki,2,’•**, 

k  ,,,  k.  k~  _,  . ..,  k  _  must  be  remembered,  i.e.,  1  +  4m  quantities, 
xuj  <*  J.  5  $  &  j  3  m-i  •  j 

In  this  case,  therefore,  if  N  =  12,  m  =  3;  so  a  system  of  three  equations 
could  be  handled  by  the  ENIAC  if  the  Runge-Kutta  fourth  order  method  is 
used  and  13  registers  are  available.  This  method  is  therefore  superior 
to  the  previously  described  methods  for  machines  of  limited  storage  space. 

It  will  be  observed  that  the  functions  f^  are  used  four  times  in  going 

from  to  xr+1  as  compared  to  once  by  the  Moulton  and  similar  methods. 

Usually  the  formation  of  these  functions  is  the  most  difficult  part  of 
the  problem  and  thus  the  Runge-Kutta  method  is  in  most  cases  more  tedious 
for  hand  computation.  However,  for  machines  such  as  the  ENIAC  where  com¬ 
puting  time  may  be  a  small  fraction  of  total  time  spent  on  a  problem, 
and  where  it  is  nearly  as  easy  to  instruct  it  to  form  f^  four  times  as 
one,  this  is  no  drawback. 

Another  method  which  is  somewhat  similar  to  the  Runge-Kutta  method 
in  adaptability  to  high  speed  machines  but  superior  in  regard  to  space, 
is  described  below.  We  have  not  found  it  in  the  literature,  (although 
it  may  have  been  known  to  Gregory,  Euler,  Newton,  «or  Lagrange).  There¬ 
fore  we  shall  derive  it  in  detail.  In  this  development  we  shall  always 
assume  that  i  rvns  from  1  to  m.  Let  us  assume  that  the  functions  f. 

L  — 

are  analytic;  then  y^  are  also  analytic.  Let  us  denote  y^(xn  +  j)  by 

_  dy.  h  _ 

V  +  7>  *  7i'<  8tc-  Er'“dine  Jijn. 


about  x 


h  +  2 


we  have: 


;7  Milne,  and  Bateman,  Op.  cit.,  pn.  77-80 
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(4.9)4 


^i,n+l  =  V  V  7  -  4’  r +  yi"'  vs  *  °(h5> 

yi,n  *  yi  -  Vs  +  yi"  {r  -  V  *5  +  yi""  4  +  °<h5) 
n,nd  ■  yi  +  *i  5  +  yi”'  5-  ♦  yi”  ZB  * 0(h4) 


7i,n 


IV.- 

where  0(hn)/hn  is  a  quantity  which  approaches  a  finite  limit  as  h 
approaches  0. 

Adding  and  subtracting  the  first  and  second  pairs  of  equations: 

,3 


-  -  h  -  -  "H  V.3  / 

=  Ji'  -  ?!*  3  +  y±”  %  -  y i  +  Js  +  °(i>  )> 


(4.io)-4 


a)  yi,Ml  -  yi,n  =  yih  +  yi'  7Z  +  0(h5) 

b)  yi,n+l  +  51,1.  =  2yi  +  5?  r  +  0(h4) 

c)  yi,n+l  -  yi,n  =  yi  h  +  °^3) 


^  yi, n+1  +  yi, n  "  2yi  +  V  r  +  0(h4)  ’ 

Assume  now  that  y.  and  x  are  known  and  that  y.  -  is  known  to 

onier  3;  i.e.,  1<”  n  1-n+1 

yi,n+l  =  yi,n+lJ  +  0(h2  1) 

where  j  =  0,  1,  2,  or  3.  Then  the  following  sequence  of  computations 
will  71.14  7ijn+1J+1 

fa)  yi,n  =  fi<yl,n’  y2,n-  — -  ym,n  xn>  equa-  <4-« 


b)  yi,n+l,3“  fi(yl,n+l,3’"-’ym,n+l,j’  Vh)+  0(h'’cq^>  (4-1) 

o)  y*1>ah  =  yi>n+1,j  -  yi,n  +  “O-h  .1“*-  (4.10)c 

2yi,j  “  yi,n+l,j  +  yi,n  “  ^,3-2  b  ^4  +  °^b  kqua.  (4.10)b 

8>  yi,3  *  <yi,d-  y2,3*  ~*,J-  *n  +  «  +  (4.1) 

!)  %1-Z  h2/4  =  yi,mi,3+  yi,n  -  2yi,3  +  °‘hd+1>  <4-10)d 

=  *•»  +  ^  +  ^‘2  h3/24tOCh3+\ua.  (4.11)a 
By  repeating  this  sequence  of  operation  4  times  y^  n+1  will  be  obtained 
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to  fourth  order.  If  a  higher  order  approximation  is  desired,  it  is 
necessary  tc  adjoin  more  points,  or  differentiate  the  differential  equa¬ 
tions  to  have  more  equations .  lince  fourth  order  is  a  convenient 
approximation,  we  shall  not  discuss  the  other  possibilities. 

let  us  compare  this  method  now  with  the  previously  discussed  ones . 
Here,  as  in  the  ftunge-Kutta  method,  there  is  no  difference  between  the 
step  from  x^  to  x^  and  the  step  from  x^  to  x^^;  nor  is  there  any  diffi¬ 
culty  in  changing  the  interval.  As  for  the  tax  on  memory,  the  moment 
of  greatest  storage  requirement  occurs  when  x  ,  y.  ,  y.  ..  .  and 

n  2  j 

y.  .  must  all  be  held  in  order  oo  form  yi"  .  _  h  /4  (y-  ,n  •  may  be 
substituted  for  y.  ,  .  ,  as  scon  as  formed).  Thus  the  inequality 

/  j— X  ^ 

1  +  3m  =  :!  must  be  satisfied.  For  the  Hi  a  AC,  therefore,  k  =  4;  i.e., 
it  is  possible  to  approximate  to  fourth  order  a  system  of  four  equations. 


The  price  for  being  able  to  solve  4  instead  of  3  equations  is  the 
formation  of  f^  20  times  instead  of  4  2s  in  the  Runge-Kutta  method;  just 

as  the  price  for  being  able  to  solve  3  instead  of  2  equations  is  the 
formation  of  f  4  times  instead  of  one. 


If  the  process  described  above  is  stopped  at  the  second  iteration, 
it  amounts  to  the  Ileun  method: 

(4-U)  5  [jri(yl,n’"’3n,n’xn)'l'ft<5l,n','hfl{yX,n’  "’y,n,n’::n}’ 


Tliis  method  yields  only  a  2nd  order  approximation  to  the  solution  but 
has  the  advantage  that  only  1  v  2m  registers  are  required  for  dead  stor¬ 
age,  so  that  a  system  ox  6  equations  may  be  ’nandled  on  the  RiHAC.  This 
is  actually  th.e  method  which  lias  often  been  used  in  the  past  for  the 
I'i'jAC.  It  is  clear,  however,  that  if  four  or  less  equations  are  in¬ 
volved,  or  if  eight  or  less  are  involved  and  only  five  significant  fig¬ 
ures  are  to  be  carried  through  the  computation,  then  it  is  highly  desira¬ 
ble  to  use  the  fourth  order  approximation  described  above;  for  larger 
3 tops  may  be  used  for  a  given  truncation  error,  and  there  will  be  corre¬ 
spondingly  less  round-off  error  (only  the  round-off  error  of  the  last 
iteration  counts).  The  following  table  summarizes  approximate  esti¬ 
mates  for  various  methods  of  solving  ordinary  differential  equations 
with  the  emphasis  on  computing  machines  of  small  memory  capacity. 
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Lie  t  hod 

Moulton 
Adams - 
Bashforth 
Stephens 
Milne 

Runge- 
,  Kutta 

Clip- 

pinger 

Heun 

Analytic 

Continuation 

Succes- 
ive  approxi¬ 
mation 

Hand 

computation 

Excellent 

Good 

Poor 

Pair 

Not  usually 
practical 

Occasionally 

good 

Machine 

computation 

Poor 

Good 

Good 

Good 

Not  usually 
applicable 

Not  adapted 
for  machines 
of  small 
memory 

Different 
procedure 
at  start 

Yes 

No 

No 

No 

No 

Interpolates 
to  reduce 
interval 

Yes 

No 

No 

No 

No 

Number  of 

registers  in  ,  .  _ 
internal  mem¬ 
ory  required 
for  storage 
for  system  of 
m  equations  and 

4th  order 
approximation 

1  +  4m 

l+3m 

(l+2m,  but 
gives  2nd 
order  approxi¬ 
mation,  not 

4th) 

Number  of 
equations 

0, 

3 

4 

(6,  but 

that  can  be  2nd  order) 

solved  to 
4th  order 
with  ENLAC 

Higher  Order  Approximations  to  the  Solution  of  Systems  of  Hyperbolic 
Partial'  bifferentlal  Equations  ‘  ■ 

Just  as  it  was  possible  to  approximate  the  solu^  on  of  a  system  of 
ordinary  differential  equations  by  a  polynomial  of  arbitrary  order  over 
a  given  interval,  so  it  is  possible  to  do  the  same  for  a  system  of  par¬ 
tial  differential  equations .  We  shall  not  attempt  to  give  any  general 
theory;  however,  certain  general  observations  may  be  made.  Suppose  we 
have  a  system  of  m-2  quasi-linear  hyperbolic  first  order  partial  dif¬ 
ferential  equations  in  two  independent  variables  u^  art  Uj,  and  m-2 

depandent  variables  u^,  uffi.  Suppose  there  are  m-2  families  of 

real  characteristic  curves,  only  two  families  being  distinct;  then 
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characteristic  variables  <X  and  may  be  introduced  so  that  the  system 
takes  the  canonical  form: 


(4.13) 


rm 

au. 

2 

— i-  =  b. 

3a  J 

i=l 

m 

3u. 

2 

i=l 

aij 

—2.  -  b. 
d/3  3 

j  =1,  2,  . ..,  & 


j  m  . 


If  now  the  a.  .  and  b.  are  analytic  functions  of  il.  ,  ....  u  and  the 
ij  J  J  J-  m 

u's  are  all  known  at  a  set  of  n  points  P^,  ?£>  ate.,  in  the  rt  ,/ff 

plane,  the  problem  of  finding  the  u's  at  a  point  P  nearby  may  be  attacked 
as  follows:  Assuming  that  the  u's  are  analytic,  vie  may  expand  them  at 
P  ,  P„,  ...,  p  ,  etc.,  about  some  convenient  point  P.  These  mn  equations 
i  *  n  _  v+l 

plus  the  m  differential  equations  at  P  plus  the  (2  -  2)m  differential 

equations  at  P  obtained  by  k-fold  differentiation  yield  ^k+l  +  n_i)m 

algebraic  or  transcendental  equations  for  u,  a^/dot  ,  &\/d/3  , 

2—  2— 

3  ,  3  ,etc.  They  may  in  general  be  solved 


Figure  4.1 


k+1 

for  (2  +n-l)m  of  these  derivatives  in  terms  of  the  known  u's  at  P^, 

P„, ...P  .  If  these  then  be  substituted  in  the  Taylor  expansions  for  the 
i  n 

u^  at  P  about  p,  a  polynomial  approximation  of  some  order  j  is  obtained. 
Since  there  are  m(1+^,..+j)  =  mj(j+l)/2 

derivatives  of  u's  at  P  of  order  j  or  less,  j  would  in  general  be  the 
largest  interger  less  than  or  equal  to 


J  =  -i  + 


>/ 9+8(2k+1 


+n-2) 


38 


For  example  if  the  differential  equations  are  not  differentiated  and 
two  points  are  used,  J  =  1  and  we  are  assured  of  being  able  to  obtain 
a  first  order  approximation  to  the  u's  at  P. 

However,  it  may  well  be  that  by  choosing  the  points  P^,  P2,  etc., 

properly  the  equations  obtained  in  {u,  3 ,  etc.,  contain  fewer 

independent  quantities  to  be  eliminated  than  mj(j+l)/2.  For  example, 
using  four  points  P^,  P2,  P y  P  arranged  at  the  corners  of  a  rectangle 

with  sides  parallel  to  the  d\  ,  /3  axes  and  with  P  chosen  at  the  center, 
the  expressions 

2-  2- 

0  B  u.  0  3  u. 

5  +  (Aa)2 - i-  +  (A  fi)2 - i 

3a  a/3 

occur  in  the  expansion  of  u^  at  each  of  the  four  points ;  thus  the  num¬ 
ber  of  quantities  of  order  two  or  less  at  P  to  be  eliminated  is  reduced 
from  6  to  4  and  so,  using  only  the  original  differential  equations  a 
second  order  approximation  to  the  u's  at  P  may  be  found.  Until  the 
reduction  in  number  of  independent  quantities  was  noted,  one  would  have 
predicted  that  there  were  only  enough  equations  to  provide  a  first  order 
approximation  to  the  u's  at  P. 

In  the  process  of  extending  the  functions  u^  from  points  P^,  ...,  ?n 

to  point-  P  we  have  used  (with  reservations)  expansions  of  u.  at  Pn 

—  k+1  1 

about  P  and  the  differential  equations  and  their  2  -2  sets  of  derived 

__ 

equations  at  P.  (2  -l)mn  additional  equations  in  the  derivatives  of 

u.  at  P  may  be  obtained  by  expanding  the  derivatives  of  u.  up  to  order 

1  —  k+1  1 

k  about  P  and  substituting  in  the  (2  -l)m  differential  equations  at 

each  point  P^,  P2,  ...,  Pn>  Except  for  checking,  these  equations,  how¬ 
ever,  are  usually  of  academic  interest  only,  since  they  have  essentially 
been  used  already  in  the  earlier  determination  of  the  u's  at  P^,  P2,  etc. 

The  assumption  that  the  u's  are  analytic  must  be  considered  for  each 

problem.  It  is  well  known  that  the  solution  of  a  system  of  hyperbolic 

differential  equations  with  analytic  coefficients  geed  n^t  be  analytic. 

(E.g.,  the  two-dimensional  wave  equation  $  u/gx  =  3  u/^y  is 

equivalent  to  the  system  p„  =  q  ,  p  =  q  ,  where  p  =  u  ,  q  =  u_ .  The 

x  y  y  jc  x 

general  solution  is  p  =  f'  (x+y)  +  g'(x-y),  q  =  f '  (x+y)-g'  (x--. ,  where 
f(ot  )  and  g (j3)  are  any  functions  with  continuous  second  derivatives, 
and  f'(ct)  =df/dc^  ,  g'(^)  =  dg/d^.)  In  fact,  the  characteristics 
may  be  defined  as  curves  along  which  ‘discontinuities  of  derivatives  of 
some  order  may  occur  even  though  the  solutions  y.  are  themselves  con¬ 
tinuous.  As  an  example  for  the  aerodynamicist,  Consider  the  flow  over 
a  body  of  revolution  with  a  contour  having  discontinuous  slopes  as  in 
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Figure  4.2.  Consider  any  curve  ABC  intersecting  at  a  non-zero  angle  the 
^  characteristics  which  bound  the  expansion  regions.  The  velocity  com¬ 
ponents  and  the  pressure  and  density  of  the  air  all  have  discontinuous 

derivatives  at  P, ,  P,.,  P_,  P,  with  respect  to  arc  length  along  this  curve. 

~  ^  / 


Figure  4.2 

On  the  other  hand,  if  functions  u^  have  analytic  initial  values  given 

along  an  analytic  curve  which  is  not  a  characteristic,  then  the  u’s 
will  be  analytic  in  the  region  of  detemLnacy.  More  general  theorems  of 
this  nature  are  available. 


Figure  4.3 

Let  us  return  now  to  the  general  process  described  in  section  3, 
and  rediscuss  it  in  the  spirit  of  the  above  general  remarks.  As  in 
section  3>  we  assume  that  we  are  to  determine  a  solution  of  equations 


a)  (K-R)yflt  -  Lxa  =  0  or  ya  =  X  xa 


b)  Hy.  -  (K-It)x_  =  0  or  (ojA  =  x 


j  r  *  -i 

\o)  Hu0e+(K-R)v0{  +  xa^-^r~+  (BRytfg)/2j  =0 


4.4  R.Courai'  andD.  Hilbert,  Ifethoden  der  Mathematischen  Physik,  vol.  II, 
Chapter  5. 

F.  Frankl  and  P.  Aleksieva,  Op.  cit.,  p.  793. 

H.  Lewy,  Op.  cit.,  p.  179* 
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r~  2 

d)  (K-R)uj  +  Lvg  +  y^  je~  -  (BRy^g)/2  J 

e)  ZY  =  “  y^Av  z  =  yeku 
with  boundary  conditions 

2 


=  0 


(4.14)4  .  Cu2  ~  ql)J _ 

cont'd  e  ~  TTT  TT7~  T 
(u2-b)(V  q-X/V6) 


du2 

35" 


v22(b-u2)  =  (q1  -  U2)2(U2  -  d) 

^  =  (ql  -  u2)/v2 

y  =  F(x) 
v  =  uF‘  (x) 

l,-o 

We  assume  that  x,  y,  u,  v,  z  are  kncmat  Jt  and  u  in  the  Ch  ,  (3  plane 
and  are  to  be  found  at  b,  the  intersection  of  characteristics  through 
Jt  and  u.  Let  us  translate  the  oi  ,  axes  to  P,  the  midpoint  of  A 
and_u.  Let  us  expand  x  at  ji  ,  b,  and  u  about  P,  indicating  quantities 
at  P  bv  bars  and  quantities  at  b  by  no  subscript,  A  oc  and  A  /3  by  oi 
and  ft  respectively,  derivatives  with  respect  to  ot  by  subscript  Oi  , 
etc. : 


x'-x+oix^  +  ftxg  +  i(oc 


2- 


x«o  +  ' 


^/3/3)+0i^Oc/3+  °(h3) 

(4.15)4  x^  =  x  -<xxa  +Ax/3  +  l-(«2xaa+A2xg;fl)-c^x0^  +  0(h3) 


U  =  x  +«xa-/3^+  lia2xa0(+/32x 


-****/£ 


where  h  is  the  larger  of  oc  and  .  Adding  and  subtracting  the  third 
equation  to  and  from  the  second,  we  obtain 


2- 


.2- 


(4.16)< 


x^  +  xu  =  2x  +  «%a+X3  x^-  20 fix^  0(h  ) 

vt  +  yu  =  2y  +  «  %***%-  W  °<h3> 

Uj  +  uu  =  2u  +  a  2U^+  /32^-  2a^*3+  <i(h3) 

Vi  +  Vu  =  ^  ■’  a  2Va+^2^"  0(h3) 

2r  .  *2-  _  -  2o^i^+  0(h3/ 


z«  +  z  =  2z 
.*  u 


+  «  *M+^ 


145 
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xd  =  xi  "  xu  =  ‘  2oi*U  +  2&*fl  +  0(h3) 

yd  =  VJL  ~  yu  =  "  2U  y*  +  7/3  +  °(h3) 

(4.17)-  ud  =  Ug  -  uu  =  -  2«  +  2/Sug  +  0(h3) 

vd=vJt  ~  \  =  "  2aV  +  2t^A  + 

Zd  =  z£  ~  \  =  "  2<*  Soc  +  2/3  Z/5  +  0(h3^ 

Substituting  from  equations  (4*14)  into  equations  (4*17),  we  ob¬ 
tain  the  two  pairs  of  equations 

xd  =“  2ax«+  *2/*%  +  °th3) 

yd  =  -  A2«x  +  2 fly  +  0(h3) 

(4*18 )J  d  *  f  _  -  _  3 

]ud  =  -  [C*P-5)/(K-R)J  2/5^  -  2«u(Jt-A2^v4+  0(10 

(vd  =  [(wWI)/(R-5)]  2«xfi(-(-E32«uoi  +  2/Svg  +  0(h3) 

for  the  pairs  of  unknowns  2<xxa  ,  2/3  and  2«u^  ,  2/3vg  .  A3  for 
the  coefficients,  we  observe  that  from  equations  (4.16) 


f(x,y,S,v,5) 


=  f  (■ 


[X-  +x 
— 2 — 

x.  +x  y  4y 
fi _ u  _Ju 


+  0(h2),  0(h2 ),..., +  0(h2)j 

L,...,  VjL_  )  +  0(h2) 


if  f  is  analytic  (class  C^).  therefore  to  the  order  indicated  the  coef¬ 
ficients  may  be  replaced  xby  fun_xions  of  the  means  of  x,  y,  u,  v,  z 
at  A  and  u.  The  determinant  of  the  coefficients  in  each  pair  of  equa¬ 
tions  is 


1  -  15*  2a  ; 

uv  +  a  /q2  -  a2 

it  would  only  be  zero  if  the  Ifech  number  were  one  or  infinity,  and  in- 

2  2 

q  -a  '  =  p.  We  rule  out  these  three  cases.  With  this 
restriction,  equations  (4.18)  may  be  solved  with  third  older  errors: 

f2  #  =  [2>yd  -  *d]  /  L1  "  ^  ^]+  °Ch3) 

(4.20X  2/6  jrg  =  j^yd  -  A*d]  /[l  -  )w]+  0(h3) 
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(4.20)- 

(cont'd) 


=  -  £(1  -  Xo>)(Xvd  +  ud)  +  {§-T -X  W  (S+T jpd  + 
2TXxd]/  [i- Xw]2  +  0(h3) 

=  [(1  -  \  w  )(vd  +  (iud)  -  2Twyd  +  ^s+f-X  w(S-T)j- xd 

/  [l-  X  w]2  +  0(h3) 


where 


e  _  €a  V 

s  “yO^E) 


and  T  = 


ByeRg 

2HPS7  * 


These  quantities  may  now  be  substituted  into  the  expansions  for  x,  y,  u, 
v,  z;  but  the  results 


X  =  X  +  J  [2  wyd  -  (1  +  X  w)xd]  /  [1  -X  coj 
y  =  y  +  2  |ji+Xa>)yd  -  2^xd]/[i  -X  wj 
(4.21)  U  =  u  +  I  [{(1+  AS  )ud  +  2  X  Vj}  -  2  {(y-y^OS-T)-, 

X  (x~xt  )  (S+T)}]/  [  1- X  U>] 

v  =  v  +  i[{2wud  -i(l+Xo»)vd}~  2  {(x-x£  )(S+f)-t5(y-yu). 

(S-T)}]/  [l-Xio] 
z  =  z  +  yA  £-(x-x)v  +(y-y)u  ]  , 

where  x  =  -|(x^  4  xj,  y  =  j(y£  +  yj,  u  =  \(u£  +  uu),  v=%-(v£  -tvj, 

f*  O  1 

are  correct  only  to  first  order  since  £[«.' X<xa  +/®  x/3/9  +  2  ex 

etc.,  have  not  been  evaluated.  These  expressions  are  the  same  as  those 
obtained  in  the  last  section.  However,  it  is  possible  to  obtain  expres¬ 
sions  which  are  correct  to  second  order  by  the  simple  expedient  of  using 
the  values  of  u^  at  a  in  Figure  4.3*  In  fact  corresponding  to  equations 

(4.16)  are  equations 


'x— x&  =  2ocxa  +  2/3xg  +  0(h3) 

y-y  =2*ya  +  2/3y«  +  o(h3) 
(4.22W  a  01  #  . 

u-ca  =  20(11^  +  2/3^0  +  0(lr) 

v-va  =  2  *v^  +  2/3^  +  0(h3) 

^z-za  =  2*^  +  2/3  z^  +  0(h3) 
which  with  equations  (4.20)  yield  equations 
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X  =  xa  +  [2^yd  -(1+ A  £)xdJ/[I- A  £>]  +  o(h3) 
y  -  ya  +  [  (i+Xw)yd  -  2Axd]/£i-X  u>]  +  o(h3) 
u  =  \  +  [{(l+Xw)ud+2Xvd}  -  2  {  (y-yu)(S-T)  - 

(4*23)  j  \(x-xx  )(3+T)}]/[  1-X Sb]  +0(h3) 

v  =  \  +[{2Wud+(l+X^)vd}  -2  {(x-x^  )(S+T)  - 

®(yyu)Cs-f)}j/£i-X£3]  +o(h3) 
z^  =  z&  +  yK  ^-v(x-xa)  +  u(y-ya)J  +  0(h3) 


It  is  important  to  note  that  the  computations  involved  in  the  second 
order  approximations  are  almost  identical  with  those  in  the  first  order 
approximations.  Accordingly,  the  only  excuse  for  using  the  first  order 
ones  would  be  in  the  case  where  storage  memory  in  a  machine  was  too  small 
to  remember  the  values  at  the  extra  point  a. 


An  alternative  procedure  for  obtaining  a  second  order  approximation 
to  x,  y,  u,  v,  z  at  b  which  is  slightly  better  for  the  machine  of  small 
memory  may  be  found  as  follows :  one  might  expect  that  if  he  replaced 

20 by  x-x^  #  2a  by  y-y^  ,  etc.,  in  equations  (4.14)a  and 
(4.14)c  evaluating  coefficients  at  the  midpoint  of  JL  and  b,  and 

similarly  replaced  by  x-x  ,  etc.,  in  equations  (4.3A)b  and  (4.14)d 

dOC  U 

evaluating  coefficients  at  the  midpoint  of  b  and  u,  a  second  order 
approximation  would  result;  and  this  i3  true.  To  prove  it,  let  us  form 
the  expressions 


(4.24) <  _ 


y  -  y i  -  X(x-x^)  =  2  u/3  A^  x^+  o(h3) 

X  -  xu  -  <5(y  -  yu)  =  2a/3  y^  +  0(h3) 

£  (u-u^  )+(v-v^  )+C3+T)(x-xtf  )  =  -2  <X/3  [dig  na  +(3+T)^  5  1 
+  Ofh3')  J 


(u-uu)+  X(v-vu)+(S-?)(y-yu)  +(^')<x  7^1 


+  0(h3) 

and  replace  X  t  CO  ,  X  >  etc.,  using  equations 

A=  Xj+  A  u)  +  °(h  ) 

W=  |(CO«  +  )  +  0(h2) 

*  u 


(4.25)-<  (S+T)  =  \  [(Si  +  Tj  )  +  (3U  +  Tu) J  +  0(h2) 

<*5^  =  z(x-x£  )  +  0(h2) 
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<1 =  f(y-yu)  +  o(h2) 

|c<Xa=  5<  A  -Ai)  4  0(h2) 


etc. 

As  a  result  we  find  the  equations 
r 


C4.26)K 


(y-y^  )  -  2'C  A4-  )(x  -  x^ )  =  o(h3) 

(x-xu)  -  1C  W  +  u)(y  -  yu)  =  0(h3) 

(V-Vj  )  +  k<*>  +  w£)(u-u/  )  +  I  |"(3+T)+CSi  +  Tj  )|  . 

(x-x£  )  =  0(h3)  L 

(U“V  +  arX  +  Au)(v-vu)  +  I  f(S-T)+(Su-Tu)J  • 

Cy-yu)  =  o(h3) 

These  may  be  solved  for  x,  y,  u,  v,  z. 

r 

x  =  [(v^6)+wu)yul  +  4<X+X^)xz}]  / 

[1  -  £(> +  X/)(" +<*>u  )] 

y  =  [{yi  “2'CA+X^  )xy}  +|(  X+Xy )  {ru-J( W  +  ^u)yu}]  / 

[i  -  i(X  +  Xi)(  wu)] 

(4‘27)  =  [{V^(  A  +  Xu)vu-Lij(y-yu)}  -  |( X  +  X  u)  {l(w  +«/ )  up  + 

▼i  -^(x-Xj  )}]  /  [l  -|(X+ X  u)(w+toi)] 

V  =  [{!(<*> +  6)/^  +  V£  -  I^Cx-x^  )}  -£(  W+  wi)  {uu+i(  X  +  Xu)vu 

-^(y-yu)}]/l_14<X+Xu)(  oj+uj£  )] 

z  =  \Jz(zx^jl  )]t {2(yA)+(yA)u  +Cja)|  }  {-(2v+vu  ^ i  )* 
(^"V1  i  )+(2u+uu+uJ)(2y-yu-y ^  )}] 

where  1^  =  |  [CS+T)+(S^  +T^  )j  and  1^4  [(S-T)+(Su-Tu)J  ,* 

but  these  formulae  must  be  used  by  iteration  to  yield  second  order  approxi¬ 
mations  since  X,  a>  ,  P,  Q,  Ayu,  Ayv  must  themselves  be  known  to  first 
order. 


We  turn  our  attention  now  to  a  third  order  approximation  to  x,  y,  u, 
v,  z  at  a  point  not  on  a  boundary.  As  we  have  suggested,  we  may  increase 
the  order  of  approximation  by  adjoining  more  points  or  by  differentiating 
the  differential  equations  or  both.  Because  of  the  fact  that  we  want  our 
methods  to  apply  when  there  are  characteristics  which  are  lines  of  dis¬ 
continuity  of  derivatives  of  x,  y,  u,  v,  z,  we  prefer  not  to  adjoin  more 
points.  The  reason  is  that  using  only  the  points  J?  ,  u,  a  and  b  on  the 
vertices  of  a  rectangle  we  can  always  manage  to  have  isolated  character- 
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istics  of  discontinuity  be  part  of  our  network  and  thus  avoid  having 
them  pass  through  any  of  the  rectangles .  It  will  then  be  permissible 
always  to  assume  that  the  dependent  functions  are  analytic  inside  the 
rectangle . 


Following  the  general  remarks  above,  we  observe  that  using  x,  y,  u, 
v,  z  as  given  at  three  arbitrary  points  and  expanded  about  a  fourth_F 
we  would  have  fifty  unknowns  (x,  xa  ,  x^  ,  x*«*  ,  x^  ,  ,x««to<, 

k+1 

etc!)  to  determine  from  the  15  +  5  (2  -1)  expansions  and  differential 

equations  obtained  by  differentiating  k  times.  Thus  k  would  have  to  be 
2.  The  problem  of  differentiating  the  differential  equations  twice, 
obtaining  35  differential  equations,  and  solving  these  35  with  15 
expansions  for  the  fifty  unknowns  would  appear  to  be  quite  formidable. 
Actually,  because  of  the  symietry  of  the  four  points,  a,  b,  jt  ,  u,  with 
respect  to  P,  the  midpoint,  the  number  of  unknowns  is  so  reduced  that  we 
shall  only  need  to  use  live  of  the  ten  possible  fi  rst  derived  equations . 


To  be  pfeeise*  if  we  expand  any  function  p  such  as  x,  y,  u,  v,  z,\, 
fed  j  ebe»j  hi  dj  b,  2  ,  u,  taking  the  origin  at  P, 

$  =  p  +c*p*  +p?s  +  *■  s(<*  3w 

3^’  )+  ^(3 01  +/8  Ppfifl  ) 

^  =  p-*^  V?  %  +  ifo2?**  +/s2p^  ^  3Pc««w  + 

(4.28)  <  ) 

pu  =  p  +«.p^  -/SPe  +i(°<.2p0^+/fl2p^  3p *«<*.+ 

Pa  =  P  -*P* -/S  P^  +-^2p^  )+<yJP^  -  5^3W  + 

we  may  form  the  cambihat.ionsi 


(4.29) 

P  +  Pa- 

■Pjt  "  Pu 

=  +  O(h^)  | 

(4.30) 

P/+Pu 

=  2p  +  0(n  ) 

(4.31) 

*P«  = 

(p  +  ptt  - 

Pa  “  P i  )/4  +0(h3)  , 

(4.32) 

ll 

ip5 

(p  +  p,  - 

P..  -  Pu)/4  +  0(h3)  , 

ci  U 

(4.33) 

P-Pa  = 

=  2«p,  + : 

+  0(h3)  ,  and 

(4.34) 

pj"pu 

+  2/fl  P5  +  0(h3)e 

Equation  (4.29)  will  yield  third  order  approximations  to  x,  y,  u,  v, 
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) 


z  provided  we  can  evaluate*^  x*^  ,  o(/Q  7^^  >  etc.  This  suggests 

that  we  differentiate  equations  (4.14)a),  c),  e)  and  (4.14 )b),  d)  with 
respect  to  /3  and  oc.  respectively  and  solve  for  4*^5^-  ,  4 o^/By^g, 
etc:  P 

4oC/S  5 =[<2  2/S  2*5*+  2*03*  2$yfi]/  1  -A^J 

4oC/d  7*0  =[\  2oCo^ 2/370  +  2<X  xj/[l  ~A  £>] 

J  4oC/?  =  U<5+f  )4of^  Tfyg  -(S-T  )4oQ5  7W/a  +S  2  ^6^2 

(4*35)]  2**2^]  /[l_-j$ 

[A 2/3  (3+T),  2*5*  -  2<*(5-T)*  2/3,  y^j  [l->  ®] 

=[uS  (S-T)4«C£yW/tf  -(3+T)4oC/S  x*^  +co  2cc>2  #  - 

-  2/3<ofl2ocuot]  /  [l  -A &]+ 

[0  2«(S-f)at  2/9^  -2£  (S+T)^  2«  y  /  [l  ->*] 

4 «/fli0(/9=  yA(-v  4*/Sxa<<s  +  u  4-C/S^  )  +  2/3 (yAu)^  2*^- 

^  2/0  (yAv)^  2«(xo(  • 

In  order  that°*/0  x*^  shall  be  correct  to  third  order  it  is  necessary 

that  A  ,u>  ,  P,  Q,  yAu,  yAv  be  known  to  first  order,  and*  S*p  > 

etc.,  and<*Aw  ,  &  A ^  ,  o<  Pw  ,  etc.  be  known  to  second  order.  The 

first  set  may  be  found  by  using  equation  (4*32).  The  second  set  may  be 
found  by  solving  equations  (4*18). 

2ocxoc  =[w  (y£  -  ya)  -  (Xjg  -  xjl  /  [l  -Aco] 

2&fi  =  [(3^-7u)  ‘*&] 

(4.36)  2*.uec  =  -[jcujj  -  uu)  +  >  ty-  vu)+3(yjg-yu)  -  T(y-ya)J  / 

[i-Aw] 

2^/9  =  [w(ux  -  uu)  +  (v^  -  vu)  +  3(x£  -xu)  -  T(x  -  xa)]  / 

[i  -A*] 

To  find  the  third  set  we  use  equations  (4.23)  and  (4*33)  and  (4.34) • 

Second  and  third  order  boundary  processes  may  also  be  devised  al¬ 
though  we  do  not  have  a  third  order  process  yet  which  is  elegant  enough 
to  .include  in  this  report. 

.To  obtain  a  second  order  contour  process,  for  example,  let  us  call  a 
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and  b  two  points  on  the  contour  in  the  oc /8  plane,  the  point  of 
intersection  of  vertical  characteristic  through  a  and  horizontal  charac¬ 
teristic  through  b,  P  and  P,  the  midpoints  of  ,  b  and  a,  b  respectively. 
We  assume  that  x,  y,~u,  v,  z  are  known  at  Ji  and  a  and  are  desired  at 
b.  We  take  the  origin  at  P,  expand  x,  y,  u,  v,  z  at  J.  and  b  about  P  , 
and  at  a  and  b  about  P,  and  form  the  differences  and  sums. 


i.  P*  *b 

P. 


Figure  4.4 


(4.37) 

P 

“pi 

=  2a  p 

a  +  o(h3), 

(4.38) 

P 

+  pi: 

=  2p  +a 

2p  +  0(h4), 

£  oca  v  '  ’ 

(4.39) 

P 

"  Pa  ' 

"  2aPa 

+  2/3^  +  0(h3), 

(4.40) 

P 

+  Pa  : 

=  2p  +a 

denoting  x. 

y» 

n,  v. 

z,  \  , 

etc.,  generically  by  p. 

+  o(h4). 


may 


Now  if  the  map  of  the  contour  is  /Q  ■-  Ta, 
be  written 


(4.41) 


dy  ^  4  r*g 

K  xa  + 


and 


the  boundary  conditions 


(4.42;  v  =  u  f(x,y)  . 

Multiplying  the  second  member  of  equation  (4*41)  above  and  below  by 
2  /3  /T  or  2  OC  we  get 


(4.43)  2aya  +  2/fly^  =  f(x,y)  [}  <X  x^+  2 /<?x^]  . 

Equations  (4.42)  and  (4.43)  must  hold  in  particular  at  ?  and  then  we 
find  the  equations 


(4.44) 


ry  -  ya  -  ?  (x  -  xa)  =  0(h3) 
r-Yj  ~  A(x-  x,)  =  0(h3) 
a>(u  -  u^)  +  (v  -  v^)  +  (S+T)(x  -  x^)  =  0(h3) 

v  -  u  f(x,y)  =  0(h3) 

V- 
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which  may  be  solved  for  second  order  approximation  to  x,  y,  u,  v,  z 
provided  X.  £  ?.  However,  would  only  equal  ?  if  the  Mach  angle  were 
zero,  i.e.  the  Ikch  number  infinite.  Ruling  this  case  out,  X  will  not 
equal  f  if  the  grid  size  is  small  enough  because  of  the  continuity  of  \  , 
These  equations  are  iterated  to  get  the  second  order;  i.e.,  X  ,  f,  etc., 
are  first  taken  equal  to  X^  »  fa,  etc.,  for  a  first  order  approximation, 

then  X  ,  f,  etc.,  are  found  from  equations  (4«38)  and  (4*40)  and  used 
in  equations  (4.44)  for  the  second  order  approximation. 


Alternatively,  in  regions  where  extrapolation  is  permissible  ^ 
J,  etc.,  may  be  found  to  first  order  by  extrapolation. 


Section  5«  Error  Study 
Introduction 


Numerical  solution  of  differential  equations  involves  making  approxi  - 
mations  which  introduce  errors.  An  estimate  of  these  errors  is  of  great 
importance  in  the  analysis  of  the  computations. 

The  first  and  thus  far  the  only  type  of  body  studied  by  the  methods 
described  in  the  preceding  sections  is  the  cone -cylinder.  A  report  on 
these  computations  will  be  published  later.  Before  they  were  carried 
out  on  the  ENIAC,  a  study  of  the  errors  involved  in  the  numerical  pro¬ 
cesses  was  made.  It  is  expected  that  the  behavior  of  the  errors  in  the 
cone-cylinder  problem  typifies  the  behavior  of  errors  in  the  calculation 
of  flows  about  arbitrary  pointed  bodies  of  revolution.  The  method  of 
analysis  and  the  differential  equations  are  the  same,  only  the  boundary 
conditions  are  specialized. 

SHIAC  Computations,  Empirical  3tudy 

To  investigate  the  effect  of  grid  size  and  order  of  approximation  on 
the  computations,  the  ^.on  in  the  expansion  region  for  a  particular  cone- 
cylinder  and  Mach  number  was  computed  by  the  EIJIAC.  The  case  studied  was 
one  for  which  9  =  20°  and  =  2.12966.  The  1st  order,  2nd  order  itera¬ 

tive,  and  2nd  Jrder-3  ooint  methods  were  used  for  grid  sizes  varying  from 
h=l  to  h=l/4C. 

The  expansion  region  in  the  physical  plane  (ABC)  is  shown  in  Figure 
5.1 ;  and  in  the  characteristic  plane  (AA'BC),  in  Figure  5.2.  (h  =  1 
corresponds  to  the  grid  oize  which  yields  x,  y,  u,  v  at  C  in  one  step 
(see  Figure  5.1).  V/hen  h  =  four  steps  are  required  to  obtain  x, 
y.  u,  v  at  C;  v/hen  h'=  l/3>  nine  steps  are  required;  etc.)  Yfe  assume 
that  the  input  data  (x,  y,  u,  v  along  the  characteristics  cX  =  0  and 
ft  -  0)  contain  no  error. 


A 


.Figure  5.1  Figure  5.2 

It  is  possible,  modifying  the  method  used  by  Frankl  and  Aleksieva"’ 

5 .2 

or  Goursat^*  ,  to  obtain  limits  on  the  size  of  the  error  in  the  values 
of  x,  y,  u,  v  at  any  point  in  the  flow  field.  These  limits  are  of  con¬ 
siderable  theoretical  interest.  However,  in  order  to  obtain  general  re¬ 
sults  it  is  necessary  to  make  rather  strong  assumptions.  The  limits  so 
obtained  are  much  larger  than  necessary,  and  give  only  a  poor  idea  of 
the  behavior  of  the  error  as  a  function  of  grid  size  or  order  of  approxi- 

nation.  The  work  of  Richardson' J  suggests  that  a  much  more  exact  study 
of  the  error  is  possible  for  each  specific  flow  problem.  A  natural  pro¬ 
cedure  is  to  solve  the  problem  for  more  than  one  grid  size,  h;  and  then 
at  points  common  to  several  grids  to  fit.  x,  y,  u,  v  to  some  reasonable 
function  of  h.  Assuming  that  these  functions  are  valid  approximations 
for  all  sufficiently  small  values  of  h,  it  is  possible  to  extrapolate 
for  the  limits  of  x,  y,  u,  v  as  h  approaches  zero. 

The  above-mentioned  procedure  was  applied  to  the  SHIAC  computations 
of  the  cone-cylinder  expansion  region.  The  computed  functions  were 
plotted  against  the  grid  size  h  at  the  points  (1,§),  (i,l),  and 

(1,1)  in  Figures  5*3,  5*4,  5*5,  and  5.6  respectively  for  the  1st,  2nd 
order  iterative,  and  2nd  order-3  point  methods.  The  curves  for  the  2nd 
order  computations  were  drawn  with  zero  slope  at  h  =0. 

The  curves  drawn  at  these  four  representative  points  should  indi¬ 
cate  the  general  behavior  of  the  error  in  the  whole  region.  The  graph¬ 
ical  extrapolations  to  h  =  o  are  shown  in  Table  5*1.  Also  included 
for  comparison  are  the  computed  values  for  h  =  1/32  by  the  2nd  order 
iterative  method. 


5.1  Frankl  and  P.  Aleksieva,  Op.  Cit.  ref.  1.9 

5.2  Goursat,  Cours  d»  Analyse,  Par.  386 

5.3  L.  F.  Richardson,  "ftie  deferred  approach  to  the  limit,  part  I  - 
single  lattice®,  Phil.  Trans.,  vol.  226,  1927,  p.  299 
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Table  5.1 

(d=l,  ,$=  2)  x  y  u  v 


1st  ord. 

mwm 

2nd  ord.  iter. 

i. 833 88 

.705713 

.723392 

-.019765 

2id  ord. -3  pt. 

1.83396 

.70565 

.723400 

-.01978 

2nd  ord.iter.,h  =  1/32 

1.83391 

.705693 

.723393 

-.019777 

(d  =  £-,  i) 

X 

y 

u 

V 

^-■wKSseasum 

2nd  ord.  iter. 

1.707045 

.764781 

.662496 

-.081612 

2nd  ord. -3  pt. 

1.70736 

.764710 

.66248 

.08161 

2nd  ord.  iter.,h=l/32 

1.707056 

.764761 

.662483 

.081609 

». 

ii 

II 

X 

y 

u 

V 

Mouam 

j masssiWiMt?: 

|  1 1 

2nd  ord.  iter. 

3.20525 

1.80904 

0.669896 

.043910 

2nd  ord. -3  pt. 

3.20550 

1.80891 

0.66978 

.04390 

2nd  ord.  iter.,h=l/32 

3.20531 

1.80899 

0.669783 

.043918 

(ct=  1,j6=  1) 

X 

y 

u 

V 

iim«^ 

r  MW.m  " 

2nd  ord.  iter. 

3.71593 

1.54547 

.709311 

2nd  ord. -3  pt. 

3.7162 

1.54535 

.70932 

-.02630 

2nd  ord.  iter.,h=l/32 

3.71610 

1.54535 

.709290 

-.026319 

The  agreement  of  the  graphical  extrapolations  indicates  that  we  may 
be  assured  of  the  accuracy  of  x,  y,  u  to  within  3  in  the  fourth  figure 
and  of  v  to  within  1  in  the  third  figure  in  all  the  values  listed  in 
Table  5.1.  Ignoring  the  first  order  extrapolation,  we  find  even  closer 
agreement, which  indicates  that  the  second  order  iterative  computation 
for  the  grid  size  h=l/32  is  reliable  to  within  6  in  the  fifth  figure 
for  x,  y,  u  and  to  within  2  in  the  4th  figure  for  v. 

The  graphs  show  that  the  2nd  order-3  point  calculations  follow  two 
different  patterns  in  relation  to  grid  size.  The  points  lie  on  either 
one  of  two  curves,  depending  on  whether  l/h  is  odd  or  even  for  (<*=1, 
f&  =1),  and  on  whether  l/2h  is  odd  or  even  for  the  other  three  points. 

To  obtain  a  numerical  extrapolation  to  h  =  0  one  can  fit  the  data 
by  the  least  squares  method  to  some  reasonable  function  of  h.  This  wa3 
done  for  the  computations  at  (&  =  1,  =  1);  the  functions  employed 

are  f(h)  =  a  +  bh  +  ch^  for  the  1st  order  values,  and  f(h)  =  a  +  bh  + 

2  q 

ch  +  dir  for  the  2nd  order  iterative  and  2nd  order-3  point  values. 

Table  5 .II  shows  the  results. 
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2nd  0R0-3PT 


SIZE  ON  CALCULA' 


Table  5. II 


1st  order 

x  =  3.71556035  -  .340695693  h  +  .046337492  h2  5 

y  =  1.54562504  +  .181446572  h  -  .00554870548  hZ 

u  =  .70931271  +  .02430564  h  -  .00444836206?h 

v  =-.026314228  -  .048620732  h  +  .011581095  h 


2nd  order-iterative 

x  =  3.715892542  +  .002385437  h  +  .134694176  h? 

y  =  1.545510530  -  .001357115  h  -  .117729376  h% 

u  =  .709308737  +  .000075481  h  -  .020277732  h  „ 

v  =-.02630000  -  .00012539676  h  -  .00125911740  li 


-  .002492252  hi 
+  .025024049  K 
+  .013444669  h£ 

-  .009568309  h3 


2nd  order-3  point 

x  =  3.716199722  -  .000608705  h  -  .002697709  h2  +  .008292560  h3 
y  =  1.545365028  +  .0003562452  h  +  .0040463651  hZ  +  .0051117669  h3 
u  =  .709310875  -  .00005547079  h  +  .0096804297  h*  -  .0077314754  h3 
v  =  -.0263339776  +  .0001389884  h  -  .05993978635  h  +  .0205477726  h3 


The  constant  terms  In  the  polynomials  are  the  extrapolations  to 
h  =  0.  The  range  of  grids  used  was  h  =  ^  to  h  =  1/40.  Only  even  grids7, 
were  used  for  the  2nd  order-3  point  data  fit.  The  above  polynomials 
all  agree  with  the  given  data  at  least  to  within  5  in  the  6th  figure 
for  x,  y,  uj  and  to  within  1  in  the  3rd  figure  for  v.  (The  v  agree¬ 
ment  for  the  first  order  computations  holds  only  for  h  £  1/20.)  Com¬ 
parison  shows  that  this  numerical  extrapolation  to  h  =  0  agrees  with 
the  graphical  extrapolation  to  within  3  in  the  5th  figure  for  x,  y,  u; 
and  to  within  1  in  the  3rd. figure  for  v. 

These  results  suggest  that  the  error  can  be  moderately  well  repre¬ 
sented  by  simpler  functions j  namely,  by  bh,  for  the  1st  order  method, 
and  by  bh^  for  the  2nd  order  methods.  The  least  squares  fits  of  the 
data  to  these  functions  are  given  in  Table  5 .III. 

Table  5. in 

1st  order 

x  =  3.715007365  -  .3284527749  h 

y  =  1.545690071  +  .1799814167  h 

u  =  .7093652534  +  .02312716466  h 

v  =  -.02644997335  -  .04280403863  h 

2nd  ord.  iter. 

x  =  3.716010825  +  .1426135703  h2 

y  =  1.545421374  -  .1164789567  li  „ 

u  =  .7093000262  -  .01673578393  k 

v  =  -.02631802511  -  .0023052374  h 


Table  5. Ill  (cont'd) 


2nd  ord,-3  pt. 

x  =  3.716162678 
y  =  1.545376882 
u  =  .7093153981 
v  =  -.0263463811 


-  .002?846l7l6L,h2 
+  .0065855559  hi 

+  .0075518560  K 

-  .0543126060  h 


These  curves  differ  from  those  in  Table  II  at  most  by  5  in  the 
5th  figure  for  x,  y,  uj  and  tjy  1  in  the  3rd  figure  for  v. 


Hound-0 ff  Errors 


Although  the  outputs  of  the  ENIAC  in  this  problem  contain  ten  fig¬ 
ures,  they  have  less  than  ten  significant  figures,  for  several  reasons. 
In  order  to  make  allowance  for  the  wide  ranges  of  sane  of  the  quantities 
encountered,  certain  numbers  such  as  x  and  y  were  purposely  shifted  to 
the  right  on  the  accumulators .  While  this  procedure  insured  that  no 
numbers  would  exceed  the  capacity  of  the  machine  in  the  extreme  cases, 
it  meant  a  loss  of  one  or  two  significant  figures.  Furthermore  ENIAC 
multiplication  and  division  are  only  correct  to  nine  places.  Thus  local 
computations  are  affected  by  round-off  errors  in  the  sixth  or  seventh 
significant  figure  of  x  and  y,  and  the  eighth  significant  figure  of  u. 

As  for  v,  since  its  magnitude  is  small  (it  can  change  sign)  it  may  have 
between  zero  and  eight  significant  figures  locally  correct.  When  it 
has  none,  however,  it  does  not  affect  the  accuracy  of  the  other  quanti¬ 
ties. 


These  local  errors,  of  the  round-off  variety,  are  in  addition  to 
the  errors  due  to  the  replacement  of  derivatives  by  difference  quotients. 
It  is  the  principal  aim  of  our  error  study  to  determine  empirically  the 
nature  of  these  latter  truncation  errors  ■> 

Because  of  the  above-mentioned  round-off  errors  we  cannot  hope  by 
extrapolation  to  zero  grid  size  to  obtain  more  than  six  significant  fig¬ 
ures. 


Theoretical  Study  of  the  Truncation  Error 


It  is  natural  to  expect  that  the  local  errors  in  computing  the  val¬ 
ues  of  x,  y,  u,  v,  at  b  (assuming  correct  values  given  at  1,  a,  and  u) 
are  cf  order  j  (i.e.  the  error  is  a  series  in  h  starting  with  terms  of 
order  j).  Then  the  total  . 

errors  at  B,  made  in  com-  P 

puting  x,  y,  u,  v  step-by 
step  from  boundaries  CAD  q 

can  be  expected  to  be  of 

£.  A 


u 


Figure  5.7 


D 
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the  form: 

(5.1)  x-x.  =  E  =  a  hn  +  E  h134^  c  hn+2+  .... 

Jj  X  X  X  X 

y-yL  =  Ey=  ay  hn  +  Ey  h3341  +  ....,  etc. 
n  =  Di 

where  x^  is  the  limit  of  x  as  h  approaches  zero,  and  the  coefficients 
a  ,  E  ,  c  ,  etc.  are  functions  of  a  and  A  independent  of  h,  hut  dif- 
fering  for  x,  y,  u,  v. 

If  this  is  true  and  if  this  series  converges  rapidly  enough  for 
the  values  of  h  for  which  the  flow  problem  has  been  solved,  we  may 
neglect  all  but  the  first  one  or  two  terms  and  solve  for  n,  a,  b  and 
x^  (or  y^  etc.).  For  example,  if  t  is  used  generically  for  x,  y,  u 

or  v,  and  if  t^,  t and  t^  are  the  values  of  t  when  h  is  h, ,  h and 

hy  we  should  have  the  equations 

(5.2)  tx  =  tL  +  l  hf 

tn  =  t.  +  a  h.n 

*3  =  h  +  *  h3n 
for  t^,  a,  and  n  to  satisfy. 

Eliminating  a  and  tT ,  n  must  satisfy  the  equation 

(5.3)  (t2-t3)  -  hj”)  =  (t1-t2)(lu,n  -  hjn). 


In  particular  if  ^1  _  ^2  h„ 

T  T  3 

(5.4)  r  t  -  t  “i 

n=LioeY^-j /  ^ 2 

n  =  3.32194  log1(J 

Having  found  n,  a  and  t^  are  given  by 


(5.5) 


a  = 


\  ~  *2 
u  n  v,  n 

ki  -h2 
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run 


(5.6)  -  a 


v  ^ 

Howe/ver,  for  some  cases  T - j— —  may  be  negative,  and  therefore 

t2"  *3 

the  data  cannot  be  fitted  for  three  values  of  h  by  a  function  of  the  form 

t  =  t.  +  a  hn» 

L 

This  means  that  the  term  5  h  cannot  be  neglected.  If  the  approx- 
mat  ion 

(5.7)  t  =  tL  +  ahn  +  B  hn+1, 

is  used,  the  four  values  of  h  and  t  available  are  substituted 

f*l  ■  *1  *  S  hi"  +  E  hlI1+1 

u  .  r  u  n  ,  c  ,  r 


(5.8)  -p  "  *L  +  *  h2  +  S  *2 

t?  =  tL  +  a  h3n  +  S  h3n+1 

t  =  tT  +  a  h,n  +  E  hp1 
\^4  L  4  4 

and  eliminated,  the  equations 

fV‘2  =  A  !  =  1  +  E  <V+1  ‘  ’'a'*1) 

(5.9) -  t2-t3  »  A  2  >*  (hjH-hj11)  +  E  (h2n*':L  -  hj1*1) 

x.  j.  _  a  _  r  /u  n  ,  m  .  c  n+1  .  n+l\ 


result . 


t3-t4  =  A3  =  5  (h3n-h4n)  +  S  (b^x  -  h4™) 


Therefore,  n  must  satisfy  the  equation 


(5.10; 


.  n  ,  n  .  n+1  ,n  +  1 

\  "i 

.  n  .  n  ,  n+1  .  n  +  1 

h2  -  h3  *2  "  ^ 

u  n  ,n  .  nfl  .n+1 


P  h3n"h4n 


u  nr.i  V,  J 

**3  "h4 


In  particular,  if  h2  h3  h,  ,  equation  (5.10)  becomes 

s “  =  r  =  r=  4 


(5.11) 
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or  (5.12)  A  x  -  3a2  2n  +  2A3  (2n)2  =  0. 
The  solutions  of  this  equation  are 


r  1  r  ( 3A2  -  \J 9A~  -  8A,  A  ,  A 

'•13)  n  =  [lAog.2j  [log  V  ■  ^  - - - - - -  j}> 


provided  9 A  2  ,  and 


;A2t\/ 


9A2  - 8AxA3 


Somewhat  more  generally,  if 


>0  ’ 

tu  h2  h3 

“*y  =  “*r  “  —  ■  —  h , ,  the  equation 
IT  k2  k  4 


to  be  satisfied  by  n  is 


(5.14)  A1  -  (k+1)  A  2  kn  +  k  A3(kn)2  =  0, 

r  Iff  (ML)*,  +  \A+i)2  a  22-4kA,  A  ^ 

(5.15)  n  =  [lA°S  kj  [*«(. - i-rrz- - i—  =3—3  y 

When  n  has  been  obtained,  t_ ,  a,  and  E  are  found  from  three  of  the 
linear  equations  (5 .8). 

The  value  of  equation  (5*13)  for  hand  computations  is  very  small, 
h^  must  be  made  small  enough  so  that  it  is  permissible  to  truncate 

equation  (5.1)  at  the  second  term.  Then  A3  will  be  so  small  that  it 

can  be  known  accurately  only  if  very  many  figures  are  carried. 


In  summary,  if  assumption  (5.1)  is  correct,  it  is  possible  to 
determine  n,  a,  E,  etc.  if  the  solution  is  carried  out  at  enough  grid 
sizes.  If  the  computations  are  done  by  hand,  a  tiny  error,  whether  by 
mistake  or  round-off,  affects  the  value  of  n  (as  given  by  equations 
(5*4),  (5*13)  >  or  similar  formulae)  30  markedly  that  the  study  cannot 
be  very  valuable .  In  fact,  it  is  possible  for  a  small  variation  in 
the  fifth  figure  of  the  data  at  small  grid  sizes  te  change  the  sign 

2 

of  the  discriminant  (9A  2  -  SA-^A^)  in  equation  (5 *13)  from  positive 
to  negative  to  yield  a  complex  value  of  n. 


Our  study  was  made  on  the  ENIAC,  a  machine  which  carries  ten 
figures  and  rarely  makes  mistakes.  Even  these  ENIAC  computations,  how¬ 
ever,  cannot  all  be  relied  upon  to  calculate  satisfactory  values  of  n, 
for  reasons  discussed  above  and  in  the  section  on  round-off  errors. 
Since  the  computations  of  u  were  found  to  have  the  most  significant 
figures,  we  calculated  n  with  them.  The  values  of  n,  found  by  equa¬ 
tion  (5.13)  from  the  u  data  at  h  =  1/4,  1/8,  1/16,  1/32,  are  listed  in 


61 


Table  5. IV 


Table  5. IV 


(a  =i,  &  =1/2) 

ex=l/2,^=l/2) 

(<*  =1/2 /=1)  (d  =1,0=1) 

1st  ord  .99 

1.0 

1.0  .97 

2nd  ord  -  iter.  2.1 

2.1 

2.1  2.1 

2nd  ord  -  3  pt.  1.9 

2.0 

2.0  2.0 

From  Table  5 .IV  we  see  that  the  order  of  the  gross  error  is  approxi¬ 
mately  one  for  the  1st  order  computations  and  approximately  two  for  the 
2nd  order  computations . 

If  assumption  (5.1)  is  correct,  as  the  data  seem  to  indicate,  we 
may  make  the  following  observations.  Let  us  consider  methods  of  two 
different  orders  for  computing  t,  calling  the  results  t.  and  t„.  VTe 
can  represent  the  errors  by  1 

_  n, 

(5.16)  tj_  -  tL  =  a;L  h  1  (1  4^) 

- 

t2  -  t^  =  h  (1  where  C-^  s.nd62  go  to  zero  as 

h  goes  to  zero  and  n_>n_  .  Then  for  any  ru  there  exists  an  h  =  H,  such 
that  1 

n. 

a  (ru,)  h0?  <  a  (n^)  h  1 

and  h^-H.  It  is  reasonable  to  expect  that  n(j)  is  an  increasing 
function  of  the  local  order.  Therefore,  a  method  of  any  local  order 
j  gives  more  accurate  results  than  methods  of  smaller  local  order  for 
all  grid  sizes,  h,  small  enough  (  <  IT) . 

Under  the  same  assumptions,  there  would  also  exist  some  (possibly 

smaller)  grid  size  H* ,  such  that  the  method  of  local  order  j  gives  re¬ 
sults  of  specified  accuracy  with  less  total  labor  than  any  method  of 

lower  local  order.  This  is  true  as  long  as  the  specified  accuracy  isf 

as  good  as,  or  better  than,  the  accuracy  associated  with  grid  size  K  . 

Although  this  conclusion  ignores  the  effects  of  round-off  errors,  it 
is  probably  correct  even  with  round-off  errors,  provided  enough  signifi¬ 
cant  figures  are  carried.  The  computation  with  higher  order  local 
approximation  has  fewer  round-off  errors  (since  it  is  carried  out  at 
larger  grid  size)  if  the  specified  accuracy  is  hign  enough. 

On  the  other  hand,  for  grid  sizes  larger  than  U*  there  ’.rill  be 
lower  order  methods  which  give  more  accurate  results  for  a  given  amount 
of  labor.  For  this  reason  and  for  reasons  of  accuracy  in  extrapolating 
to  zero  grid  size,  the  second  order  method  was  found  the  best  for  hand 
and  machine  computations. 
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Having  found  the  2nd  order  iterative  method  most  feasible  for 
computation  tvith  the  ENIAC  we  wish  to  investigate  the  accuracy  obtain¬ 
able  by  the  process  of  computing  at  large  grid  sizes  and  extrapolating 
to  h  =  0*  This  method  involves  considerably  less  labor  than  computa¬ 
tion  at  very  small  grid  sizes. 


Letting  t  be  a  generic  symbol  for  x,  y,  u,  and  v,  we  write  for  t 
the  seccnd  order  function  of  h: 

t(h)  -  tL  =  ah2  , 

where  t  is  the  desired  extrapolated  value.  If  we  use  two  grid  sizes 
h^  and  such  that  h^  =  1/2  h^,  then 

(5.17)  tL  =  tChg)  +  1/3  [t(h2)  -  tC^)]  . 


Using  =  1/L6  and  =  1/32,  the  smallest  available  grid  sizes 

applicable  to  equation  (5*17),  we  calculate  t  throughout  the  expansion 

region,  and  we  consider  it  to  be  the  "correct"  value  of  t.  Then  we  form 
the  quantity  A  ^  =  £t^  -  t(h)  ^  /t^  throughout  the  expansion  region. 

A  h  is  the  relative  error  in  the  computation  at  the  point  (cX  ,  ) 

resulting  from  the  use  of  the  finite  grid  size  h. 

We  also  calculate  tj/r  =  t(l/8)  +  1/3  £t(l/8)  -  t(l/4) J  and 

form  A  *  =  (t^  -  t^)/t^.  It  is  A  *  that  we  wish  to  compare  with  A  ^ 

for  various  grid  sizes  to  see  how  tliu  accuracy  of  extrapolation  with 
large  grids  compares  with  that  of  small  grid  computations. 


Ir.  figure  5.8  we  have  the  relative  errors  in  x  plotted  for  the  grid 
sizes  l/8,  l/l6,  1/32,  and  for  the  1/4,  1/8  extrapolation.  It  is  evi¬ 
dent  that  the  l/4,  1/8  extrapolations  are  not  as  good  as  the  1/32  x  1/32 
computations  but  are  appreciably  better  than  the  l/l6  x  l/l6  computations. 
The  errors  in  y,  u,  and  v  behave  in  an  identical  manner  with  those  of  x. 

2 

The  amount  of  labor  required  for  grid  size  h  is  proportional  to  l/h  . 
Taking  the  extrapolations  as  equivalent  to  computations  with  h  =  l/28, 

2  2  2 

the  ratio  of  the  work  required  is  about  (4+8  )/(28)  ,  approximately 
1  Ain 


Similar  conclusions  can  be  drawn  from  the  1st  order  and  2nd  order- 
3  point  method  computations. 


R.  F.  Clippinger 
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